FillNueSandbox_module.cc
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////
2 /// \brief Module to calculate and store sandbox variables in testing but
3 /// added to the CAFs.
4 /// \author rtoner@fnal.gov, bckhouse@fnal.gov
5 //////////////////////////////////////////////////////////////////////////////
6 #include <map>
7 #include <vector>
8 
9 #include "TVector3.h"
10 
14 #include "fhiclcpp/ParameterSet.h"
15 
16 #include "Calibrator/Calibrator.h"
17 #include "GeometryObjects/Geo.h"
19 #include "Geometry/LiveGeometry.h"
20 #include "MCCheater/BackTracker.h"
21 #include "NueSandbox/NueSandObj.h"
22 #include "RecoBase/CellHit.h"
23 #include "RecoBase/RecoHit.h"
24 #include "RecoBase/Cluster.h"
25 #include "RecoBase/Track.h"
26 #include "Utilities/AssociationUtil.h"
27 #include <TParticlePDG.h>
28 #include <TDatabasePDG.h>
29 
30 namespace nuesand
31 {
33  {
34  public:
35  explicit FillNueSandbox(const fhicl::ParameterSet& pset);
36  virtual ~FillNueSandbox();
37 
38  void produce(art::Event& evt);
39  void beginRun(art::Run& run);
40  double GetMincell(const rb::Track& track);
41  int GetMincellWall();
42  double GetTrackAngle(const art::Ptr<rb::Track> track);
43  void GetRegression(const float sumz, const float sumz2, const float sumr, const float sumzr, const float totweight, float& slope, float& inter);
44  double TrackD(const float pz0,const float pr,const float slope,const float inter);
45 
46  /// Fraction of plane pairs where the "track" seems to change angle, ie
47  /// have a different number of hits/plane. For a real track this will be
48  /// low, for showers it should be higher.
49  ///
50  /// loose=true allows 1 cell per plane without counting it in the score
51  double FracAngChanges(const rb::Cluster* slice, bool loose) const;
52 
53  /// Fraction of planes having the modal number of hits/plane. For a
54  /// straight track this should be high, for a shower lower.
55  double FracModalHits(const rb::Cluster* slice) const;
56 
57  /// Fraction of planes where the cells aren't all adjacent. For a track
58  /// this should be low, for a shower or two-pronged event, high.
59  ///
60  /// loose=true allows single cell gaps and still call a plane contiguous
61  double FracNonContiguous(const rb::Cluster* slice, bool loose) const;
62 
63  /// Maximum number of contiguous unoccupied planes. Should be low for
64  /// normal events, high for misslicings (and possibly pizeros).
65  int MaxGap(const rb::Cluster* slice) const;
66 
67  /// Number of hits in most-populated plane of \a view
68  int MaxHits(const rb::Cluster* slice, geo::View_t view) const;
69 
70  /// Fraction of energy of the mother neutrino contained inside the detector
71  /// In the NC case, the energy of the daughter neutrino is ignored
72  double GetECF(const art::Ptr<simb::MCTruth> mctruth) const;
73 
74  /// Calculate the energy in plane
75  double GetPlaneEnergy(const art::PtrVector<rb::CellHit>& chits,
76  const art::Ptr<rb::Prong>& prong,
77  const int & iPlane);
78 
79  /// Compute dE/dx of a prong, averaged over all planes
80  /// Does not count 0's for planes with no hits in them.
81  double GetdEdx(const art::Ptr<rb::Prong> prong);
82 
83 
84  protected:
86  const std::vector<const rb::CellHit*>& allhits,
87  const art::Ptr<rb::Cluster> slice) const;
88 
95  float fMIPlow;
96  float fMIPhigh;
97 
99 
100  //Cluster to contain MIP hits
102 
103  private:
104  int mincellwall=0;
105 
106  };
107 
108  //----------------------------------------------------------------------
110  {
111  double fwdcell, bakcell;
112 
116 
117  TVector3 start = track.Start();
118  TVector3 end = track.Stop();
119  TVector3 dirfor = (end-start).Unit();
120  TVector3 dirbak = (start-end).Unit();
121  double point[3], dir[3];
122  int wall;
123 
124  end.GetXYZ(point); // starting at track end point projected forwards along track dir = fwd
125  dirfor.GetXYZ(dir);
126  fwdcell = livegeom->ProjectedLiveCellsToEdge(point,dir);
127  wall = livegeom->ProjectedWall(point,dir);
128 
129  start.GetXYZ(point); // starting at track start point projected backwards along track dir = bak
130  dirbak.GetXYZ(dir);
131  bakcell = livegeom->ProjectedLiveCellsToEdge(point,dir);
132 
133  double mincell = std::min(bakcell,fwdcell);
134 
135  if(fwdcell < bakcell){ // this wall is assuming the Instrumented Det Ends (downstream); if there are dropped DCMs (or other non-rectangular solid geometry) it might be wrong
136  mincellwall = wall;
137  }
138  else{
139  mincellwall = livegeom->ProjectedWall(point,dir);
140  }
141 
142  return mincell;
143 
144  }
145 
146  //----------------------------------------------------------------------
148  {
149  return mincellwall;
150  }
151 
152  //----------------------------------------------------------------------
154  {
156 
157  TVector3 trackdir = track->Dir();
158  TVector3 beamdir = geom->NuMIBeamDirection();
159  float ret = trackdir*beamdir;
160 
161  return ret;
162 
163  }
164 
165  //----------------------------------------------------------------------
166 
167  void FillNueSandbox::GetRegression(const float sumz, const float sumz2, const float sumr, const float sumzr, const float totweight, float& slope, float& inter)
168  {
169 
170  slope = 0;
171  inter = 0;
172 
173  if (totweight>0){
174 
175  if (sumz2 == (1/totweight)*sumz*sumz){
176  slope = 1.0e4;
177  } else {
178  slope = (sumzr - (1/totweight)*sumz*sumr)/(sumz2 - (1/totweight)*sumz*sumz);
179  }
180  inter = sumr/totweight - slope*sumz/totweight;
181  }
182 
183  }
184 
185  //----------------------------------------------------------------------
186 
187  double FillNueSandbox::TrackD(const float pz0,const float pr,const float slope,const float inter)
188  {
189 
190  double dtrk = 0;
191 
192  double pz = pz0;
193  if (pz0<0) pz=0;
194 
195  dtrk = (pz + slope*(pr - inter))/TMath::Sqrt(1+slope*slope);
196 
197  if (dtrk<0) dtrk=0;
198 
199  return dtrk;
200 
201  }
202 
203  //----------------------------------------------------------------------
205  bool loose) const
206  {
207  std::map<int, int> hitsPerPlane[2];
208 
209  const unsigned int cellMax = slice->NCell();
210  for(unsigned int cellIdx = 0; cellIdx < cellMax; ++cellIdx){
211  const rb::CellHit& chit = *slice->Cell(cellIdx);
212  ++hitsPerPlane[chit.View()][chit.Plane()];
213  }
214 
215  int nChanges = 0;
216  for(geo::View_t view: {geo::kX, geo::kY}){
217  int hitsPrev = -1;
218  for(auto it: hitsPerPlane[view]){
219  const int hits = it.second;
220  if(hitsPrev > 0){
221  if(loose){
222  if(abs(hits-hitsPrev) > 1) ++nChanges;
223  }
224  else{
225  if(hits != hitsPrev) ++nChanges;
226  }
227  }
228  hitsPrev = hits;
229  }
230  }
231 
232  return nChanges/double(hitsPerPlane[0].size()+hitsPerPlane[1].size()-2);
233  }
234 
235  //----------------------------------------------------------------------
236  double FillNueSandbox::FracModalHits(const rb::Cluster* slice) const
237  {
238  std::map<int, int> hitsPerPlane[2];
239 
240  const unsigned int cellMax = slice->NCell();
241  for(unsigned int cellIdx = 0; cellIdx < cellMax; ++cellIdx){
242  const rb::CellHit& chit = *slice->Cell(cellIdx);
243  ++hitsPerPlane[chit.View()][chit.Plane()];
244  }
245 
246  double ret = 0;
247  for(geo::View_t view: {geo::kX, geo::kY}){
248  int maxCells = 0;
249  for(auto it: hitsPerPlane[view]) maxCells = std::max(it.second, maxCells);
250  int numPlanes = 0;
251  for(auto it: hitsPerPlane[view]){
252  if(it.second == maxCells) ++numPlanes;
253  }
254  ret += (.5*numPlanes)/hitsPerPlane[view].size();
255  }
256 
257  return ret;
258  }
259 
260  //----------------------------------------------------------------------
262  bool loose) const
263  {
264  std::map<int, std::vector<int> > hitsInPlane;
265 
266  const unsigned int cellMax = slice->NCell();
267  for(unsigned int cellIdx = 0; cellIdx < cellMax; ++cellIdx){
268  const rb::CellHit& chit = *slice->Cell(cellIdx);
269  hitsInPlane[chit.Plane()].push_back(chit.Cell());
270  }
271 
272  int numNonContig = 0;
273 
274  const int thresh = loose ? 2 : 1;
275 
276  for(auto planeIt: hitsInPlane){
277  std::vector<int>& hits = planeIt.second;
278  std::sort(hits.begin(), hits.end());
279 
280  int lastCell = -1;
281  for(int cell: hits){
282  if(cell-lastCell > thresh && lastCell >= 0){
283  ++numNonContig;
284  break;
285  }
286  lastCell = cell;
287  }
288  }
289 
290  return numNonContig/double(hitsInPlane.size());
291  }
292 
293  //----------------------------------------------------------------------
294  /// Maximum number of contiguous unoccupied planes. Should be low for normal
295  /// events, high for misslicings (and possibly pizeros).
296  int FillNueSandbox::MaxGap(const rb::Cluster* slice) const
297  {
298  std::set<int> planes;
299 
300  const unsigned int cellMax = slice->NCell();
301  for(unsigned int cellIdx = 0; cellIdx < cellMax; ++cellIdx){
302  planes.insert(slice->Cell(cellIdx)->Plane());
303  }
304 
305  int prevPlane = -1;
306  int maxGap = 0;
307  for(int plane: planes){
308  if(prevPlane >= 0) maxGap = std::max(maxGap, plane-prevPlane-1);
309  prevPlane = plane;
310  }
311 
312  return maxGap;
313  }
314 
315  //----------------------------------------------------------------------
317  {
318  std::map<int, int> hits_per_plane;
319  for(unsigned int i = 0; i < slice->NCell(view); ++i){
320  ++hits_per_plane[slice->Cell(view, i)->Plane()];
321  }
322  int ret = 0;
323  for(auto it: hits_per_plane) ret = std::max(it.second, ret);
324  return ret;
325  }
326 
327  //----------------------------------------------------------------------
329  fSlicerLabel (pset.get< std::string >("SlicerLabel")),
330  fCosmicLabel (pset.get< std::string >("CosmicLabel")),
331  fKalmanLabel (pset.get< std::string >("KalmanLabel")),
332  fProngLabel (pset.get< std::string >("ProngLabel")),
333  fProng3DInstLabel (pset.get< std::string >("Prong3DInstLabel")),
334  generatorName (pset.get< std::string >("GeneratorInput")),
335  fMIPlow (pset.get< float >("MIPlow")),
336  fMIPhigh (pset.get< float >("MIPhigh")),
337  fSkipPionTruth (pset.get< bool >("SkipPionTruth"))
338 
339  {
340  produces<std::vector<nuesand::NueSandObj>>();
341  // Associate with the slice
342  produces<art::Assns<nuesand::NueSandObj, rb::Cluster>>();
343 
344  produces< std::vector<rb::Track> >();
345  // Associate Tracks with the Slice they came from
346  produces< art::Assns<rb::Track, rb::Cluster> >();
347 
348  }
349 
350 
351  //----------------------------------------------------------------------
353  {
354  }
355 
356  //-------------------------------------------------------------------
358  {
359  }
360 
361  //-------------------------------------------------------------------
363  {
364 
368 
369  std::unique_ptr<std::vector<nuesand::NueSandObj>> sandcol(new std::vector<nuesand::NueSandObj>);
370  std::unique_ptr<art::Assns<nuesand::NueSandObj, rb::Cluster>> assns(new art::Assns<nuesand::NueSandObj, rb::Cluster>);
371  std::unique_ptr< std::vector<rb::Track> > trackcol(new std::vector<rb::Track>);
372  std::unique_ptr< art::Assns<rb::Track, rb::Cluster> > assnsslc(new art::Assns<rb::Track, rb::Cluster>);
373 
375  evt.getByLabel(fSlicerLabel, slices);
376 
377  art::FindManyP<rb::Track> cosmicAssnList(slices, evt, fCosmicLabel);
378  art::FindManyP<rb::Track> kalmanAssnList(slices, evt, fKalmanLabel);
379  art::FindManyP<rb::Prong> prongAssnList(slices, evt,
382 
383  const unsigned int sliceMax = slices->size();
384 
385  // Use property that every hit is in exactly one slice
386  std::vector<const rb::CellHit*> allhits;
387  for(const rb::Cluster& slice: *slices){
388  for(unsigned int i = 0; i < slice.NCell(); ++i){
389  allhits.push_back(slice.Cell(i).get());
390  }
391  }
392 
393  std::map<int,std::vector<int> > plane2cell;
394 
395  for(unsigned int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
396 
397  //Clear out map of planes->cells
398  plane2cell.clear();
399 
400  //Ptr to cluster
401  const art::Ptr<rb::Cluster> slice(slices, sliceIdx);
402 
403  //Continue if noise
404  if(slice->IsNoise()) continue;
405 
406  //Get slice mean and fid (to define sandbox)
407  const TVector3 mean = slice->MeanXYZ();
408  const bool inFid = geom->isInsideFiducialVolume(mean.X(),
409  mean.Y(),
410  mean.Z());
411 
412 
413  //For Calibration:
414  const double wx = (slice->NXCell() > 0) ? slice->W(slice->XCell(0).get()) : 0;
415  const double wy = (slice->NYCell() > 0) ? slice->W(slice->YCell(0).get()) : 0;
416 
417  //Count of hits and MIP hits:
418  double tothits=0;
419  double totmips=0;
420 
421  slicemip.Clear();
422 
423  //Counters for regression analysis:
424  float totweightx=0;
425  float totweighty=0;
426  int tothitsx=0;
427  int tothitsy=0;
428  float sumz_x = 0;
429  float sumz_y = 0;
430  float sumx = 0;
431  float sumy = 0;
432  float sumz2_x = 0;
433  float sumz2_y = 0;
434  float sumzx = 0;
435  float sumzy = 0;
436 
437  float pz, px, py;
438  float pecorr = 0;
439 
440  // variables used to find extent of planes in
441  // each cluster view
442  unsigned int nXplanes = 0, nYplanes = 0;
443  unsigned int PXhi = 0;
444  unsigned int PXlo = UINT_MAX;
445  unsigned int PYhi = 0;
446  unsigned int PYlo = UINT_MAX;
447 
448  //loop over hits in each view independently
449  for(int view = geo::kX; view <= geo::kY; ++view) {
450  const geo::View_t geoview = geo::View_t(view);
451 
452  for(unsigned int hitIdx = 0; hitIdx < slice->NCell(geoview); ++hitIdx) {
453 
454  const rb::CellHit* h = slice->Cell(geoview,hitIdx).get();
455  if(h->View() == geo::kX) {
456  if(h->Plane() > PXhi) PXhi = h->Plane();
457  if(h->Plane() < PXhi) PXlo = h->Plane();
458  }
459  if(h->View() == geo::kY) {
460  if(h->Plane() > PYhi) PYhi = h->Plane();
461  if(h->Plane() < PYhi) PYlo = h->Plane();
462  }
463  }
464  }
465 
466  nXplanes = PXhi - PXlo+1;
467  nYplanes = PYhi - PYlo+1;
468 
469 
470  //------------->Loop over hits in Slice!
471  for(unsigned int hitIdx = 0; hitIdx < slice->NCell(); ++hitIdx)
472  {
473  const art::Ptr<rb::CellHit>& chit = slice->Cell(hitIdx);
474  //Get calbrated reco hit:
475  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, chit->View() == geo::kX ? wx : wy));
476 
477  //Make sure things have been calibrated:
478  if(!rhit.IsCalibrated())
479  {
480  LOG_DEBUG("NueSandbox") << "Not calibrated?! "
481  << chit->Plane() << " " << chit->Cell();
482  continue;
483  }
484 
485  //Total Hits:
486  tothits += 1.0;
487  pecorr = rhit.PECorr();
488  pz = rhit.Z();
489  px = rhit.X();
490  py = rhit.Y();
491 
492  //For regression:
493  if (chit->View()==geo::kX){
494 
495  tothitsx++;
496  totweightx += pecorr;
497  sumz_x += pz*pecorr;
498  sumz2_x += pz*pz*pecorr;
499  sumx += px*pecorr;
500  sumzx += pz*px*pecorr;
501  } else {
502 
503  tothitsy++;
504  totweighty += pecorr;
505  sumz_y += pz*pecorr;
506  sumz2_y += pz*pz*pecorr;
507  sumy += py*pecorr;
508  sumzy += pz*py*pecorr;
509  }
510 
511  //Total MIP Hits:
512  if (rhit.PECorr()>fMIPlow&&rhit.PECorr()<fMIPhigh) {
513  totmips += 1.0;
514  slicemip.Add(chit); //Add to new slice!
515  }//end mip if
516 
517  //Plane and Cell Number:
518  int fNumPlane = chit->Plane();
519  int fNumCell = chit->Cell();
520 
521  //If there is no entry for this plane, add a new key and vector:
522  if ( plane2cell.find(fNumPlane)==plane2cell.end() )
523  plane2cell.insert(std::pair<int, std::vector<int> >(fNumPlane, std::vector<int>()));
524 
525  //Push back unique cell number for this Plane (no duplicates!):
526  if (std::find(plane2cell[fNumPlane].begin(),plane2cell[fNumPlane].end(),fNumCell) == plane2cell[fNumPlane].end())
527  plane2cell[fNumPlane].push_back(fNumCell);
528 
529  }//<----------END LOOP OVER HITS
530 
531  //New fraction of MIP Hits
532  //double fracmip=0;
533  //if (tothits>0) fracmip = totmips/tothits;
534 
535  //For regression:
536  float slopex = 0;
537  float interx = 0;
538  float slopey = 0;
539  float intery = 0;
540 
541  //Get regression values:
542  GetRegression(sumz_x,sumz2_x,sumx,sumzx,totweightx,slopex,interx);
543  GetRegression(sumz_y,sumz2_y,sumy,sumzy,totweighty,slopey,intery);
544 
545  float dtrk;
546  float dtrkx_min=1.0e6;
547  float dtrkx_max=0;
548  float dtrky_min=1.0e6;
549  float dtrky_max=0;
550  float dz_start_x=0;
551  float dz_stop_x=0;
552  float dz_start_y=0;
553  float dz_stop_y=0;
554 
555  float mincellavg = 1.0e4;
556 
557  if (tothitsx > 1 && tothitsy > 1){
558 
559  //------------->Loop over hits in Slice AGAIN (regression)
560  for(unsigned int hitIdx = 0; hitIdx < slice->NCell(); ++hitIdx)
561  {
562  const art::Ptr<rb::CellHit>& chit = slice->Cell(hitIdx);
563  //Get calbrated reco hit:
564  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, chit->View() == geo::kX ? wx : wy));
565 
566  //Make sure things have been calibrated:
567  if(!rhit.IsCalibrated())
568  {
569  LOG_DEBUG("NueSandbox") << "Not calibrated?! "
570  << chit->Plane() << " " << chit->Cell();
571  continue;
572  }
573 
574  dtrk = 0;
575  pz = rhit.Z();
576  px = rhit.X();
577  py = rhit.Y();
578 
579 
580  if (chit->View() == geo::kX) { //X VIEW
581 
582 
583  //Projected distance from z=0 to point
584  dtrk = TrackD(pz,px,slopex,interx);
585 
586  //is it a max?
587  if (dtrk>dtrkx_max) {
588  dtrkx_max = dtrk;
589  dz_stop_x = dtrk/TMath::Sqrt(1+slopex*slopex);
590  }
591 
592  //is it a min?
593  if (dtrk<dtrkx_min) {
594  dtrkx_min = dtrk;
595  dz_start_x = dtrk/TMath::Sqrt(1+slopex*slopex);
596  }
597  } else { //Y VIEW
598 
599  //Projected distance from z=0 to point
600  dtrk = TrackD(pz,py,slopey,intery);
601 
602  //is it a max?
603  if (dtrk>dtrky_max) {
604  dtrky_max = dtrk;
605  dz_stop_y = dtrk/TMath::Sqrt(1+slopey*slopey);
606  }
607 
608  //is it a min?
609  if (dtrk<dtrky_min) {
610  dtrky_min = dtrk;
611  dz_start_y = dtrk/TMath::Sqrt(1+slopey*slopey);
612  }
613  }
614 
615  }//<------END SECOND HIT LOOP
616 
617  float dz_start = TMath::Min(dz_start_x,dz_start_y);
618  float dz_stop = TMath::Max(dz_stop_x,dz_stop_y);
619 
620  TVector3 dtrkstart;
621  dtrkstart[0] = interx + slopex*dz_start;
622  dtrkstart[1] = intery + slopey*dz_start;
623  dtrkstart[2] = dz_start;
624 
625  TVector3 dtrkstop;
626  dtrkstop[0] = interx + slopex*dz_stop;
627  dtrkstop[1] = intery + slopey*dz_stop;
628  dtrkstop[2] = dz_stop;
629 
630  TVector3 dtrkdir;
631  dtrkdir[0] = slopex;
632  dtrkdir[1] = slopey;
633  dtrkdir[2] = 1;
634 
635  //Make very simple track:
636  rb::Track rtrack(*slice,dtrkstart,dtrkdir);
637  rtrack.AppendTrajectoryPoint(dtrkstop);
638 
639  //Mincell for regression line:
640  mincellavg = GetMincell(rtrack);
641 
642  //Add to vector and make association
643  trackcol->push_back(rtrack);
644  util::CreateAssn(*this,evt,*trackcol,slice,*assnsslc);
645 
646  }
647 
648 
649  //Min and Max of box around MIP frac hits:
650  TVector3 def(1000, 1000, 1000);
651  TVector3 mipmin = slicemip.NCell() == 0 ? def : slicemip.MinXYZ();
652  TVector3 mipmax = slicemip.NCell() == 0 ? -def : slicemip.MaxXYZ();
653 
654  unsigned int numplanesfull = 0; //Number of filled planes
655  unsigned int maxCellContig = 0; //Max number cells in a row in one plane
656  unsigned int maxCellSum = 0; //Max total of cells in one plane
657 
658  //------------>Iterate over planes:
659  std::map< int, std::vector<int> >::iterator planeIter;
660  for (planeIter = plane2cell.begin(); planeIter != plane2cell.end(); ++planeIter) {
661 
662  //Count planes:
663  numplanesfull++;
664 
665  //Grab vector of hits
666  std::vector<int> cellvect = planeIter->second;
667 
668  //Sort into correct order
669  std::sort(cellvect.begin(),cellvect.end());
670 
671  //Look for maximum cells in single plane:
672  if (cellvect.size()>maxCellSum)
673  maxCellSum = cellvect.size();
674 
675  //Counters for max cells in a row:
676  unsigned int maxrow=0;
677  unsigned int nrow=1;
678 
679  //------------>Iterate over cells:
680  for (unsigned int icell=0; icell<cellvect.size(); icell++){
681  if (icell!=0){
682  if (cellvect[icell] == (cellvect[icell-1]+1)){
683  nrow++;
684  if (nrow>maxrow) maxrow=nrow;
685  } else nrow=1;
686  }
687  }//------------->END LOOP OVER CELLS
688 
689  //Most in a row in this plane?
690  if (maxrow>maxCellContig)
691  maxCellContig = maxrow;
692 
693  }//<------------END LOOP OVER PLANES
694 
695  //Some dummy vars:
696  float tempcell = 1.0e4;
697  float trackangle = -5;
698  int maxtrackID = -1;
699  float maxtracklen = 0;
700  float mincellcos = 1.0e4;
701  int mincellwallcos = 0;
702 
703  if(cosmicAssnList.isValid()){
704  const std::vector< art::Ptr<rb::Track> > sliceTracksCosmic = cosmicAssnList.at(sliceIdx);
705 
706  //-----------> LOOP OVER COSMIC TRACKS
707  for(size_t trkIdx = 0; trkIdx < sliceTracksCosmic.size(); ++trkIdx){
708 
709  //Track:
710  const art::Ptr<rb::Track> track = sliceTracksCosmic[trkIdx];
711 
712  //Find projected cells to edge:
713  tempcell = GetMincell(*track);
714 
715  //Is this the closest to the edge?
716  if (tempcell < mincellcos){
717  mincellcos = tempcell;
718  mincellwallcos = GetMincellWall();
719  }
720 
721 
722  //Find longest cosmic track:
723  if (track->TotalLength() > maxtracklen) {
724  maxtracklen = track->TotalLength();
725  maxtrackID = trkIdx;
726  }
727  }//<------------END LOOP OVER TRACKS
728 
729  //CosAngle vs beam for the longest cosmic track (-5 otherwise)
730  if (maxtrackID>=0)
731  trackangle = GetTrackAngle(sliceTracksCosmic[maxtrackID]);
732 
733  }//<-------End isValid FindManyP
734 
735  float mincellkal = 1.0e4;
736 
737  if(kalmanAssnList.isValid()){
738  const std::vector< art::Ptr<rb::Track> > sliceTracksKalman = kalmanAssnList.at(sliceIdx);
739 
740  //Some dummy vars:
741  tempcell = 1.0e4;
742 
743 
744  //-----------> LOOP OVER KALMAN TRACKS
745  for(size_t trkIdx = 0; trkIdx < sliceTracksKalman.size(); ++trkIdx){
746 
747  //Track:
748  const art::Ptr<rb::Track> track = sliceTracksKalman[trkIdx];
749 
750  //Find projected cells to edge:
751  tempcell = GetMincell(*track);
752 
753  //Is this the closest to the edge?
754  if (tempcell < mincellkal) mincellkal = tempcell;
755 
756  }//<------------END LOOP OVER TRACKS
757 
758  }//<----End loop over Kalman isValid
759 
760  double dedx_png1 = 0, dedx_png2 = 0;
761  if(prongAssnList.isValid()){
762  const std::vector< art::Ptr<rb::Prong> > sliceProngs =
763  prongAssnList.at(sliceIdx);
764 
765  double maxlength = 0, secMaxlength = 0;
766  int maxidx = -1, secmaxidx = -1;
767 
768  for(unsigned int ipng = 0; ipng < sliceProngs.size(); ++ipng){
769  double len = sliceProngs[ipng]->TotalLength();
770 
771  if(len > maxlength ){
772  secMaxlength = maxlength;
773  maxlength = len;
774  secmaxidx = maxidx;
775  maxidx = ipng;
776  }
777  if(len > secMaxlength && len < maxlength){
778  secMaxlength = len;
779  secmaxidx = ipng;
780  }
781  }// end loop over prongs
782 
783  if(maxidx > -1 )
784  dedx_png1 = GetdEdx(sliceProngs[maxidx]);
785  if(secmaxidx > -1 )
786  dedx_png2 = GetdEdx(sliceProngs[secmaxidx]);
787  }// endif prong assn is valid
788 
789  //Make and fill the sandbox!
790  nuesand::NueSandObj sandbox(inFid, mean);
791  sandbox.fMincellCosmic = mincellcos;
792  sandbox.fMincellWallCosmic = mincellwallcos;
793  sandbox.fMincellKalman = mincellkal;
794  sandbox.fAngleCosmic = trackangle;
795  sandbox.fNMIPHits = totmips;
796  sandbox.fMIPMin = mipmin;
797  sandbox.fMIPMax = mipmax;
798  sandbox.fFilledPlanes = numplanesfull;
799  sandbox.fMaxCellsInPlane = maxCellSum;
800  sandbox.fMaxRowOfCells = maxCellContig;
801  sandbox.fMincellAvg = mincellavg;
802  sandbox.fExtentXPlanes = nXplanes;
803  sandbox.fExtentYPlanes = nYplanes;
804 
805  sandbox.fFracAngChanges = FracAngChanges(slice.get(), false);
806  sandbox.fFracAngChangesLoose = FracAngChanges(slice.get(), true);
807  sandbox.fFracModalHits = FracModalHits(slice.get());
808  sandbox.fFracNonContig = FracNonContiguous(slice.get(), false);
809  sandbox.fFracNonContigLoose = FracNonContiguous(slice.get(), true);
810  sandbox.fMaxGap = MaxGap(slice.get());
811  sandbox.fMaxHitsX = MaxHits(slice.get(), geo::kX);
812  sandbox.fMaxHitsY = MaxHits(slice.get(), geo::kY);
813 
814  sandbox.fdEdxProng1 = dedx_png1;
815  sandbox.fdEdxProng2 = dedx_png2;
816 
817  FillTruthVars(sandbox, allhits, slice);
818 
819  sandcol->push_back(sandbox);
820  util::CreateAssn(*this, evt, *sandcol, slice, *assns);
821  } // end for slice idx
822 
823  evt.put(std::move(sandcol));
824  evt.put(std::move(assns));
825  evt.put(std::move(trackcol));
826  evt.put(std::move(assnsslc));
827  }
828 
829  //-------------------------------------------------------------------
831  const std::vector<const rb::CellHit*>& allhits,
832  const art::Ptr<rb::Cluster> slice) const
833  {
835  if(!bt->HaveTruthInfo()) return;
836 
838 
839  // Ideally we wouldn't be essentially duplicating CAFMaker's truth matching
840  // here. Maybe we can have a truth sandbox that's associated with MCTruths,
841  // not slices?
842 
843  const std::vector<cheat::NeutrinoEffPur> eps = bt->SliceToNeutrinoInteractions(slice, allhits);
844  if(eps.empty()) return;
845 
846  const art::Ptr<simb::MCTruth> truth = eps[0].neutrinoInt;
847 
848  vars.ECFNu = GetECF(truth);
849 
850  const std::vector<const sim::Particle*> parts = bt->MCTruthToParticles(truth);
851 
852  // Find the biggest pions
853  int bestIdx = -1;
854  for(unsigned int primIdx = 0; primIdx < parts.size(); ++primIdx){
855  const sim::Particle* part = parts[primIdx];
856  if(part->Mother()) continue; // Not a primary
857 
858  const int pdg = part->PdgCode();
859  const double E = part->E();
860  if(pdg == 111){
861  ++vars.fNPi0;
862  if(E > vars.fPi0Energy){
863  bestIdx = primIdx;
864  vars.fPi0Energy = E;
865  }
866  }
867  if(pdg == +211){
868  ++vars.fNPiPlus;
869  if(E > vars.fPiPlusEnergy) vars.fPiPlusEnergy = E;
870  }
871  if(pdg == -211){
872  ++vars.fNPiMinus;
873  if(E > vars.fPiMinusEnergy) vars.fPiMinusEnergy = E;
874  }
875  } // end for primIdx
876 
877  // No pizeros, leave the rest zeros
878  if(bestIdx < 0) return;
879 
880  // Otherwise, take the biggest
881  const sim::Particle* pizero = parts[bestIdx];
882 
883  // TEMPORARILY DISABLING CODE BELOW
884  // TO BE REVERTED SOON
885  // The code below calls sim::Particle::Daughter() which is sometimes
886  // busted. This will be fixed in newly produced sim files, but for
887  // now, just skip the storing of this info.
888  if(fSkipPionTruth) return;
889 
890  if(pizero->NumberDaughters() != 2) return;
891 
892  const sim::ParticleNavigator& nav = bt->ParticleNavigator();
893  const sim::Particle* phot0 = nav[pizero->Daughter(0)];
894  if(phot0->PdgCode() != 22) return;
895  const sim::Particle* phot1 = nav[pizero->Daughter(1)];
896  if(phot1->PdgCode() != 22) return;
897 
898  vars.fPhot0E = phot0->E();
899  vars.fPhot1E = phot1->E();
900  vars.fOpenCos = phot0->Momentum(0).Vect().Unit().Dot(phot1->Momentum(0).Vect().Unit());
901 
902  // If the photon leaves the detector it can go a very long way
903  const TVector3 e0 = phot0->EndPosition().Vect();
904  if(geom->isInsideDetectorBigBox(e0.X(), e0.Y(), e0.Z())){
905  vars.fConvL0 = (e0-phot0->Position(0).Vect()).Mag();
906  }
907  else{
908  vars.fConvL0 = -1;
909  }
910  const TVector3 e1 = phot1->EndPosition().Vect();
911  if(geom->isInsideDetectorBigBox(e1.X(), e1.Y(), e1.Z())){
912  vars.fConvL1 = (e1-phot1->Position(0).Vect()).Mag();
913  }
914  else{
915  vars.fConvL1 = -1;
916  }
917 
918  // Make sure the biggest photon is first
919  if(vars.fPhot1E > vars.fPhot0E){
920  std::swap(vars.fPhot0E, vars.fPhot1E);
921  std::swap(vars.fConvL0, vars.fConvL1);
922  }
923  }
924 
925  //-------------------------------------------------------------------
926 
927 
929  {
930  float ECFNu(0.0);
931  float DeltaECC(0.0);
932  float DeltaENON(0.0);
933  float TrueEnergyNON(0.0);
934  float ContainmentCC(0.0);
935  float ContainmentNN(0.0);
936  float EscapedEnergyNN(0.0);
937  float EscapedEnergyCC(0.0);
939 
940  const simb::MCNeutrino& mc_neutrino = mc->GetNeutrino();
941  const simb::MCParticle& nu = mc_neutrino.Nu();
943  std::vector<const sim::Particle*> parts=bt->MCTruthToParticles(mc); ///All the primary and other "levels" daughter particles.
945 
946  ////Yep,"CUT"
947  if(abs(nu.Vx()) < geom->DetHalfWidth() && //True Vertex Limitations (Interaction Vertex !)
948  abs(nu.Vy()) < geom->DetHalfHeight() &&
949  nu.Vz() < geom->DetLength() &&
950  nu.Vz() > 0 ){
951  if(mc_neutrino.CCNC()==simb::kCC){
952  DeltaECC = 0;
953  EscapedEnergyCC = 0.0;
954  for (unsigned int j = 0; j < parts.size(); ++j){ ///Loop for All Particles
955  bool inDetStart = fabs(parts[j]->Vx(0)) <= geom->DetHalfWidth() //Check X inside detector
956  && fabs(parts[j]->Vy(0)) <= geom->DetHalfHeight() // then check Y
957  && parts[j]->Vz(0) >= 0 && parts[j]->Vz(0) <= geom->DetLength(); //then starting and End Z
958  if (!inDetStart) continue;
959  bool inDet = fabs(parts[j]->EndX()) <= geom->DetHalfWidth() ///Check End X inside detector
960  && fabs(parts[j]->EndY()) <= geom->DetHalfHeight() ///Check End Y inside detector
961  && parts[j]->EndZ() >= 0 && parts[j]->EndZ() <= geom->DetLength(); ///Check End Z inside detector
962  pnav.Add(new sim::Particle(*parts[j])); ///Keep adding each particle to the PaticleNavigator
963  if (!inDet && ///We calculate the escaped energy of all escaped particles
964  abs(parts[j]->PdgCode()) != 12 &&
965  abs(parts[j]->PdgCode()) != 14 &&
966  abs(parts[j]->PdgCode()) != 16 ){
967  for (unsigned int q = 0; q < parts[j]->NumberTrajectoryPoints(); ++q){ ///Loop for All hit Points
968  if (fabs(parts[j]->Vx(q)) > geom -> DetHalfWidth() || ///Check to see if a trahectory point is outside the detector
969  fabs(parts[j]->Vy(q)) > geom -> DetHalfHeight() ||
970  parts[j]->Vz(q) > geom->DetLength() ||
971  parts[j]-> Vz(q) < 0) {
972  DeltaECC = (parts[j]->E(q-1) - parts[j] -> Mass());
973  } ///Check to see if trajectory point is in the detector
974  } /// All Hit Point
975  EscapedEnergyCC += DeltaECC;
976  }/// Pick up non-neutrinos particles which are inside the detecter
977  }///Loop for All Particles
978  ContainmentCC = (nu.E() - EscapedEnergyCC);
979  ECFNu = (ContainmentCC/nu.E());
980  LOG_DEBUG("NueSandbox") <<"ECF(CC) is "<< ECFNu;
981  }///CC CUT
982  //////////////&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& NC &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&///////
983  else if(mc_neutrino.CCNC()==simb::kNC){
984  DeltaENON = 0;
985  EscapedEnergyNN = 0.0;
986  for (unsigned int j = 0; j < parts.size(); ++j){ ///Loop for All Particles
987  bool inDetStart = fabs(parts[j]->Vx(0)) <= geom->DetHalfWidth() //Check X inside detector
988  && fabs(parts[j]->Vy(0)) <= geom->DetHalfHeight() // then check Y
989  && parts[j]->Vz(0) >= 0 && parts[j]->Vz(0) <= geom->DetLength(); //then starting and End Z
990  if (!inDetStart) continue;
991  bool inDet = fabs(parts[j]->EndX()) <= geom->DetHalfWidth() ///Check End X inside detector
992  && fabs(parts[j]->EndY()) <= geom->DetHalfHeight() ///Check End Y inside detector
993  && parts[j]->EndZ() >= 0 && parts[j]->EndZ() <= geom->DetLength(); ///Check End Z inside detector
994  pnav.Add(new sim::Particle(*parts[j])); ///Keep adding each particle to the PaticleNavigator
995  if (!inDet && ///We calculate the escaped energy of all escaped particles except neutrinos
996  abs(parts[j]->PdgCode()) != 12 &&
997  abs(parts[j]->PdgCode()) != 14 &&
998  abs(parts[j]->PdgCode()) != 16 ){
999  for (unsigned int q = 0; q < parts[j]->NumberTrajectoryPoints(); ++q){ ///Loop for All hit Points
1000  if (fabs(parts[j]->Vx(q)) > geom -> DetHalfWidth() || ///Check to see if a trahectory point is outside the detector
1001  fabs(parts[j]->Vy(q)) > geom -> DetHalfHeight() ||
1002  parts[j]->Vz(q) > geom->DetLength() ||
1003  parts[j]-> Vz(q) < 0) {
1004  DeltaENON = (parts[j]->E(q-1) - parts[j] -> Mass());
1005  } ///Check to see if trajectory point is in the detector
1006  } /// All Hit Point
1007  EscapedEnergyNN += DeltaENON;
1008  }/// Pick up non-neutrnos particles which are inside the detecter
1009  }///Loop for All Particles
1010  int numberOfPrimaries = pnav.NumberOfPrimaries();
1011  ///Now, we calculate the ture energy by added all the primaries' energy
1012  float RealEnergyNN = 0.0;
1013  TrueEnergyNON = 0.0;
1014  for (int i = 0; i != numberOfPrimaries; ++i ){ ///Loop for the all primary particles
1015  const sim::Particle* myprimary = pnav.Primary(i);
1016  if(abs(myprimary->PdgCode()) == 12 ||
1017  abs(myprimary->PdgCode()) == 14 ||
1018  abs(myprimary->PdgCode()) == 16 ){
1019  RealEnergyNN=(myprimary->E() - myprimary->Mass());
1020  TrueEnergyNON =nu.E()-RealEnergyNN;
1021  } ///Ignoring mother neutrinos from the primaries
1022  } ///Loop for the all primary particles
1023  ContainmentNN = (TrueEnergyNON - EscapedEnergyNN);
1024  ECFNu = (ContainmentNN/nu.E());
1025  LOG_DEBUG("NueSandbox") <<"ECF(NC) is "<< ECFNu;
1026  }///NC CUT
1027  else{
1028  LOG_DEBUG("NueSandbox") <<"It is neither CC nor NC";
1029  }
1030  }
1031  return ECFNu;
1032  }// End of GetECF
1033 
1034 
1035  //-------------------------------------------------------------------
1036 
1037 
1039  const art::Ptr<rb::Prong>& prong,
1040  const int& iPlane)
1041  {
1042  double planeE = 0;
1043  int nchits = hits.size();
1044 
1045  for(int ihit = 0; ihit < nchits; ++ihit){
1046  int hitplane = hits[ihit]->Plane();
1047 
1048  if( hitplane == iPlane ){
1049  rb::RecoHit rhit = prong->RecoHit(ihit);
1050  if(rhit.IsCalibrated()){
1051  planeE = planeE + rhit.GeV();
1052  }
1053  }
1054 
1055  // hits are ordered by plane. Take advantage
1056  else if( hitplane > iPlane)
1057  return planeE;
1058 
1059  }// end loop over chits
1060 
1061  return planeE;
1062  }// End of GetPlaneEnergy
1063 
1064 
1065  //-------------------------------------------------------------------
1066 
1068  {
1070  art::PtrVector<rb::CellHit> pngHits = prong->AllCells();
1071  rb::SortByPlane(pngHits);
1072  std::set<int> planes;
1073 
1074  double dedx = 0;
1075  double dirz = fabs(prong->Dir().Z());
1076 
1077  if( pngHits.empty() )
1078  return dedx;
1079 
1080  unsigned int ncell = pngHits.size();
1081 
1082  for(unsigned int i = 0; i < ncell; ++i)
1083  planes.insert( prong->Cell(i)->Plane() );
1084 
1085  for(int iplane : planes){
1086  double celldepth = 2 * geom->Plane(iplane)->Cell(0)->HalfD();
1087  double dx = prong->TotalLength();
1088 
1089  if( dirz > 1e-6)
1090  dx = celldepth/dirz;
1091 
1092  double de = GetPlaneEnergy(pngHits, prong, iplane);
1093  dedx = dedx + de/dx;
1094 
1095  }// end loop over planes with hits
1096 
1097  return dedx/planes.size();
1098 
1099  }// End of GetdEdx
1100 
1101  //-------------------------------------------------------------------
1102 
1104 } //end namespace
const XML_Char int len
Definition: expat.h:262
double E(const int i=0) const
Definition: MCParticle.h:232
T max(const caf::Proxy< T > &a, T b)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
double FracAngChanges(const rb::Cluster *slice, bool loose) const
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Cluster.cxx:121
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
double HalfD() const
Definition: CellGeo.cxx:205
double GetTrackAngle(const art::Ptr< rb::Track > track)
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
int fMaxHitsY
How many hits does the most-occupied Y plane have?
Definition: NueSandObj.h:131
set< int >::iterator it
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:224
std::vector< NeutrinoEffPur > SliceToNeutrinoInteractions(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of neutrino interactions...
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
int Mother() const
Definition: MCParticle.h:212
geo::View_t View() const
Definition: CellHit.h:41
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
double Mass() const
Definition: MCParticle.h:238
Vertical planes which measure X.
Definition: PlaneGeo.h:28
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
float Z() const
Definition: RecoHit.h:38
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
double DetLength() const
double Mass(Resonance_t res)
resonance mass (GeV)
void abs(TH1 *hist)
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
float fConvL1
Conversion length of photon 1.
Definition: NueSandObj.h:107
void FillTruthVars(NueSandObj &vars, const std::vector< const rb::CellHit * > &allhits, const art::Ptr< rb::Cluster > slice) const
const PlaneGeo * Plane(unsigned int i) const
float fPhot1E
Energy of largest pi0&#39;s smallest photon.
Definition: NueSandObj.h:104
TVector3 MeanXYZ(rb::AveragingScheme=kDefaultScheme) const
Definition: Cluster.cxx:538
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
void produce(art::Event &evt)
int NumberDaughters() const
Definition: MCParticle.h:216
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
int ProjectedLiveCellsToEdge(TVector3 vertex, TVector3 dir)
recommended projection to use
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
double GetECF(const art::Ptr< simb::MCTruth > mctruth) const
Definition: Run.h:31
int Daughter(const int i) const
Definition: MCParticle.cxx:112
sandbox variables
Definition: NueSandObj.h:17
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
TVector3 MaxXYZ() const
Definition: Cluster.cxx:492
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
unsigned short Cell() const
Definition: CellHit.h:40
TString part[npart]
Definition: Style.C:32
double GetPlaneEnergy(const art::PtrVector< rb::CellHit > &chits, const art::Ptr< rb::Prong > &prong, const int &iPlane)
Calculate the energy in plane.
virtual double TotalLength() const
Distance along prong to reach last cell hit.
Definition: Prong.cxx:186
void hits()
Definition: readHits.C:15
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
double dx[NP][NC]
const sim::Particle * Primary(const int) const
float fFracAngChangesLoose
Definition: NueSandObj.h:125
Float_t E
Definition: plot.C:20
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
double GetdEdx(const art::Ptr< rb::Prong > prong)
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
bool isInsideDetectorBigBox(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Big Box?
Collect Geo headers and supply basic geometry functions.
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
const std::map< std::pair< std::string, std::string >, Variable > vars
const double j
Definition: BetheBloch.cxx:29
int fNPiMinus
Number of primary pi-.
Definition: NueSandObj.h:97
double GetMincell(const rb::Track &track)
double DetHalfHeight() const
TVector3 Unit() const
float fConvL0
Conversion length of photon 0.
Definition: NueSandObj.h:106
float fPiMinusEnergy
Energy of largest primary pi-.
Definition: NueSandObj.h:101
bool empty() const
Definition: PtrVector.h:336
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:308
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
int fMaxHitsX
How many hits does the most-occupied X plane have?
Definition: NueSandObj.h:130
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
double Vx(const int i=0) const
Definition: MCParticle.h:220
int ProjectedWall(TVector3 vertex, TVector3 dir)
if we assume the downstream InstrumentedDetEnds, what wall would this projected track exit from...
int MaxGap(const rb::Cluster *slice) const
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
float GeV() const
Definition: RecoHit.cxx:69
double DetHalfWidth() const
int MaxHits(const rb::Cluster *slice, geo::View_t view) const
Number of hits in most-populated plane of view.
void GetRegression(const float sumz, const float sumz2, const float sumr, const float sumzr, const float totweight, float &slope, float &inter)
T const * get() const
Definition: Ptr.h:321
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
TVector3 MinXYZ() const
Definition: Cluster.cxx:446
const int cellMax
double TrackD(const float pz0, const float pr, const float slope, const float inter)
float X() const
Definition: RecoHit.h:36
void clear()
Definition: Handle.h:236
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
void geom(int which=0)
Definition: geom.C:163
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
float Mag() const
float fPhot0E
Energy of largest pi0&#39;s largest photon.
Definition: NueSandObj.h:103
static const double celldepth
Definition: geo.cxx:46
double FracNonContiguous(const rb::Cluster *slice, bool loose) const
double Vz(const int i=0) const
Definition: MCParticle.h:222
unsigned int fExtentXPlanes
Definition: NueSandObj.h:119
int fNPi0
Number of primary pi0.
Definition: NueSandObj.h:95
float Y() const
Definition: RecoHit.h:37
float PECorr() const
Definition: RecoHit.cxx:47
FillNueSandbox(const fhicl::ParameterSet &pset)
double FracModalHits(const rb::Cluster *slice) const
T min(const caf::Proxy< T > &a, T b)
ParticleNavigator Add(const int &offset) const
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Float_t e
Definition: plot.C:35
Event generator information.
Definition: MCNeutrino.h:18
float fOpenCos
Opening angle of largest pi0&#39;s photons.
Definition: NueSandObj.h:105
virtual void Clear()
Forget about all owned cell hits.
Definition: Cluster.cxx:279
float fPi0Energy
Energy of largest primary pi0.
Definition: NueSandObj.h:99
double Vy(const int i=0) const
Definition: MCParticle.h:221
float fPiPlusEnergy
Energy of largest primary pi+.
Definition: NueSandObj.h:100
Module to calculate and store sandbox variables in testing but added to the CAFs. ...
Definition: FillSandboxes.h:6
bool isInsideFiducialVolume(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Fiducial Volume?
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).
Definition: CellHit.cxx:124
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const
unsigned int fExtentYPlanes
Definition: NueSandObj.h:120
int fNPiPlus
Number of primary pi+.
Definition: NueSandObj.h:96