FillSandbox_module.cc
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////
2 /// \brief Sandbox to make objects for numu analysis
3 /// in CAF
4 /// \author Kirk Bays
5 ////////////////////////////////////////////////////////
6 #include <string>
7 
13 #include "fhiclcpp/ParameterSet.h"
14 
15 #include "NovaDAQConventions/DAQConventions.h"
18 #include "Geometry/Geometry.h"
19 #include "Geometry/LiveGeometry.h"
21 #include "RecoBase/Cluster.h"
22 #include "RecoBase/Track.h"
23 #include "RecoBase/Vertex.h"
24 #include "RecoBase/CellHit.h"
25 #include "RecoBase/RecoHit.h"
26 #include "ReMId/ReMId.h"
27 #include "Utilities/AssociationUtil.h"
28 #include "Utilities/func/MathUtil.h"
29 #include "NumuEnergy/NumuE.h" // New
30 
31 #include "MCCheater/BackTracker.h"
32 #include "TFile.h"
33 #include "TMVA/Reader.h"
34 
35 namespace numusand {
36 
37  class FillSandbox : public art::EDProducer {
38  public:
39  explicit FillSandbox(fhicl::ParameterSet const & pset);
40  virtual ~FillSandbox();
41 
42  void produce(art::Event & evt);
43  void beginRun(art::Run& run);
44 
45  private:
46 
53 
54  };
55 }
56 
57 
58 namespace numusand
59 {
60 
61  //----------------------------------------------------------------------
63  fSlicerLabel (pset.get< std::string >("SlicerLabel")),
64  fTrackLabel (pset.get< std::string >("TrackLabel")),
65  fEnergyLabel (pset.get< std::string >("EnergyLabel")),
66  fPIDLabel (pset.get< std::string >("PIDLabel")),
67  fVtxLabel (pset.get< std::string >("VtxLabel"))
68  {
69  produces< std::vector<numusand::NumuSandObj> >();
70  produces< art::Assns<numusand::NumuSandObj, rb::Cluster> >();
71  }
72 
73 
74  //----------------------------------------------------------------------
76  {
77  }
78 
79  //-------------------------------------------------------------------
81  {
82  fNumuSandFxs = new NumuSandFxs();
83  } //end beginRun
84 
85  //-------------------------------------------------------------------
87  {
88 
90  evt.getByLabel(fSlicerLabel,slicevec);
91 
92  if(slicevec->empty()) {
93  mf::LogWarning ("No Slices")<<"No Slices in the input file";
94  return;
95  }
96 
100 
101  art::FindManyP<rb::Track> trackAssnList(slicevec, evt, fTrackLabel);
102 
103  // set up associations between slices and energy
104  art::FindManyP<numue::NumuE> eAss(slicevec,evt,fEnergyLabel);
105 
106  std::unique_ptr< std::vector<numusand::NumuSandObj> > outputObjects(new std::vector<numusand::NumuSandObj>);
107  std::unique_ptr< art::Assns<numusand::NumuSandObj, rb::Cluster> > assoc(new art::Assns<numusand::NumuSandObj, rb::Cluster>);
108 
109  for(size_t i = 0; i < slicevec->size(); ++i){
110 
111  art::Ptr<rb::Cluster> slice(slicevec,i);
112 
113  //get the total no. of cell hits
114  art::PtrVector<rb::CellHit> allslchits = slice->AllCells();
115 
116  if(slice->IsNoise()) continue;
117 
118  numusand::NumuSandObj sandbox; // initialization in constructor
119 
120  float missE;
121 
122  fNumuSandFxs->getMissingE(slice, missE);
123 
124  sandbox.SetMissingE(missE);
125 
126  sandbox.SetNDropped(livegeom->NDropouts());
127  sandbox.SetNBad(livegeom->NBadDCMBC());
128 
129  float xmin = 9999, ymin = 9999, zmin = 9999, xmax = -9999, ymax = -9999, zmax = -9999;
130  float xmin2 = 9999, ymin2 = 9999, zmin2 = 9999, xmint = 9999, ymint = 9999, zmint = 9999;
131  float xmax2 = -9999, ymax2 = -9999, zmax2 = -9999, xmaxt = -9999, ymaxt = -9999, zmaxt = -9999;
132  int ntx = 0, nty = 0;
133 
134 
135  // Temporary code. Not a great place to put these and not much time to
136  // re-order anything. Calculate vtx E within 20/40/60 cm to origin.
137  art::FindManyP<rb::Vertex> fmElastic(slicevec, evt,fVtxLabel);
138  if(fmElastic.isValid()){
139  std::vector< art::Ptr<rb::Vertex> > elastics = fmElastic.at(i);
140  if (elastics.size() == 1){
141  // Now, we can get the vtx, and loop over cells to get vars
142  TVector3 vtx = elastics[0]->GetXYZ();
143  float vtxE20 = 0;
144  float vtxE40 = 0;
145  float vtxE60 = 0;
146 
147  for (unsigned int hIdx = 0; hIdx < slice->NCell(); hIdx++){
148  rb::RecoHit rhit = slice->RecoHit(hIdx);
149  if (!rhit.IsCalibrated()) continue;
150  bool isX = slice->Cell(hIdx)->View()==geo::kX;
151  // Calcualte dist from hit to vtx.
152  // Only calculate using coordinate cell measures! (DocDB 14330)
153  double d =
154  isX ? util::sqr(rhit.X()-vtx[0]) :util::sqr(rhit.Y()-vtx[1]);
155  d += util::sqr(rhit.Z()-vtx[2]);
156  d = sqrt(d);
157  if (d < 20){
158  vtxE20 += rhit.GeV();
159  vtxE40 += rhit.GeV();
160  vtxE60 += rhit.GeV();
161  }
162  else if (d < 40){
163  vtxE40 += rhit.GeV();
164  vtxE60 += rhit.GeV();
165  }
166  else if (d < 60)
167  vtxE60 += rhit.GeV();
168  } // End loop over NCells
169  sandbox.SetVtxE20(vtxE20);
170  sandbox.SetVtxE40(vtxE40);
171  sandbox.SetVtxE60(vtxE60);
172  } // End NVtx == 1
173  } // End isValid
174 
175 
176 
177  for(unsigned int j = 0; j < slice->NCell(); ++j) {
178  const rb::CellHit* h = slice->Cell(j).get();
179  double xyz[3];
180  geom->Plane(h->Plane())->Cell(h->Cell())->GetCenter(xyz);
181  const double x = xyz[0];
182  const double y = xyz[1];
183  const double z = xyz[2];
184 
185  if(h->View()==geo::kX) {
186  if(x < xmin) { xmin2 = xmin; xmin = x;}
187  else if(x < xmin2) xmin2 = x;
188  if(z < zmin) { zmin2 = zmin; zmin = z;}
189  else if(z < zmin2) zmin2 = z;
190  if(x > xmax) { xmax2 = xmax; xmax = x;}
191  else if(x > xmax2) xmax2 = x;
192  if(z > zmax) { zmax2 = zmax; zmax = z;}
193  else if (z > zmax2) zmax2 = zmax;
194 
195  if(h->PE() > 100) { // above threshold only TODO: add tracks and recohits, use pecorr
196  ntx++;
197  if(x < xmint) xmint = x;
198  if(z < zmint) zmint = z;
199 
200  if(x > xmaxt) xmaxt = x;
201  if(z > zmaxt) zmaxt = z;
202  }
203  }
204 
205  else if(h->View()==geo::kY) {
206  if(y < ymin) { ymin2 = ymin; ymin = y;}
207  else if (y < ymin2) ymin2 = y;
208  if(z < zmin) { zmin2 = zmin; zmin = z;}
209  else if (z < zmin2) zmin2 = z;
210 
211  if(y > ymax) { ymax2 = ymax; ymax = y;}
212  else if(y > ymax2) ymax2 = y;
213  if(z > zmax) { zmax2 = zmax; zmax = z;}
214  else if(z > zmax2) zmax2 = z;
215 
216  if(h->PE() > 100) { // above threshold only TODO: add tracks and recohits, use pecorr
217  nty++;
218  if(y < ymint) ymint = y;
219  if(z < zmint) zmint = z;
220 
221  if(y > ymaxt) ymaxt = y;
222  if(z > zmaxt) zmaxt = z;
223  }
224  }
225 
226  }
227  sandbox.SetXMin2(xmin2);
228  sandbox.SetYMin2(ymin2);
229  sandbox.SetZMin2(zmin2);
230  sandbox.SetXMax2(xmax2);
231  sandbox.SetYMax2(ymax2);
232  sandbox.SetZMax2(zmax2);
233 
234  if(ntx > 0 && nty > 0) {
235  sandbox.SetXMint(xmint);
236  sandbox.SetYMint(ymint);
237  sandbox.SetZMint(zmint);
238  sandbox.SetXMaxt(xmaxt);
239  sandbox.SetYMaxt(ymaxt);
240  sandbox.SetZMaxt(zmaxt);
241  }
242 
243  else { // -5 default can be confused as a real value; +=9999 too far from real values
244  sandbox.SetXMint(-999);
245  sandbox.SetYMint(-999);
246  sandbox.SetZMint(-999);
247  sandbox.SetXMaxt(-999);
248  sandbox.SetYMaxt(-999);
249  sandbox.SetZMaxt(-999);
250  }
251 
252 
253  if(!trackAssnList.isValid()) {
254  std::cout << "NumuSandbox: No Kalman Tracks!" << std::endl;
255  outputObjects->push_back(sandbox);
256  util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
257  continue;
258  }
259 
260  const std::vector< art::Ptr<rb::Track> > sliceTracks = trackAssnList.at(i);
261 
262  art::FindOneP<remid::ReMId> remidVec(sliceTracks, evt, fPIDLabel);
263  if(!remidVec.isValid()) {
264  // std::cout << "NumuSandbox: No ReMId information!" << std::endl;
265  outputObjects->push_back(sandbox);
266  util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
267  continue;
268  }
269 
270 
271  int bestIdx = -1; // for index in list of the highest PID track
272  int pimuIdx = -5; // for index in list of a pimu decay event
273  int Idx2 = -1;
274  int idx3d = -1;
275  double bestpimuE = -999.; // keep track of how much energy due to a non-primary muon was in a track
276  // this is for if a pi->mu muon contributes light to more than once track
277  // we only want to associate with the track that saw 'most' of it.
278  // note we assume all non-primary muons are from pi->mu decay.
279  double bestval = -999., val2 = -999;
280  int nmutrks = 0;
281 
282  // ----------------------------------
283  // WORKING AREA
284 
285  float avededxtrk1 = 0.0;
286  float avededxtrk2 = 0.0;
287  float avededxtrk1Last4Cells = 0.0;
288  float avededxtrk2Last4Cells = 0.0;
289  float avededxtrk1Last6Cells = 0.0;
290  float avededxtrk2Last6Cells = 0.0;
291  float avededxtrk1Last8Cells = 0.0;
292  float avededxtrk2Last8Cells = 0.0;
293  float scattangletrk1 = 0.0;
294  float scattangletrk2 = 0.0;
295  float actvtxx = 0.0;
296  float slccalE = 0.0;
297 
298  if((sliceTracks.size() == 2)){
299  // Loop over tracks in slice
300  for(size_t j = 0; j < sliceTracks.size(); ++j){
301  const art::Ptr<rb::Track> track = sliceTracks[j];
302  art::Ptr<remid::ReMId> trackRemid = remidVec.at(j);
303 
304  if(!track->Is3D()) continue;
305 
306  if(trackRemid->Value() > bestval) { // best track so far
307  val2 = bestval;
308  Idx2 = bestIdx;
309  bestval = trackRemid->Value();
310  bestIdx = j;
311  }
312  else if(trackRemid->Value() > val2) { // not best, but better second best track
313  val2 = trackRemid->Value();
314  Idx2 = j;
315  }
316 
317  } // end loop over tracks
318 
319  // Making sure is OK
320  if((!remidVec.isValid()) ||
321  (bestval < 0.0) || (val2 < 0.0)){
322  continue;
323  }
324 
325  // Get info on muon like track
326  //if ((sliceTracks[bestval]->NCell()>0) && (bestval>0.75)){
327  // get info about this track
328  const art::Ptr<rb::Track> track1 = sliceTracks[bestIdx];
329  //
330  // Fill here with dEdx info and other variables
331  //
332  fNumuSandFxs->getAveTrackdEdx(track1,avededxtrk1);
333  sandbox.SetAvedEdxTrack1(avededxtrk1);
334 
335  fNumuSandFxs->getAveTrackdEdxLast4Cells(track1,avededxtrk1Last4Cells);
336  sandbox.SetAvedEdxTrack1Last4Cells(avededxtrk1Last4Cells);
337 
338  fNumuSandFxs->getAveTrackdEdxLast6Cells(track1,avededxtrk1Last6Cells);
339  sandbox.SetAvedEdxTrack1Last6Cells(avededxtrk1Last6Cells);
340 
341  fNumuSandFxs->getAveTrackdEdxLast8Cells(track1,avededxtrk1Last8Cells);
342  sandbox.SetAvedEdxTrack1Last8Cells(avededxtrk1Last8Cells);
343 
344  scattangletrk1 = fNumuSandFxs->getScatt(track1);
345  if(track1->TotalLength() > 0) sandbox.SetScattAngleTrack1(scattangletrk1/track1->TotalLength());
346 
347  //}
348  //else if((sliceTracks[val2]->NCell()>0)){
349  // get info about this track
350  const art::Ptr<rb::Track> track2 = sliceTracks[Idx2];
351  //
352  // Fill here with dEdx info and other variables
353  //
354  fNumuSandFxs->getAveTrackdEdx(track2,avededxtrk2);
355  sandbox.SetAvedEdxTrack2(avededxtrk2);
356 
357  fNumuSandFxs->getAveTrackdEdxLast4Cells(track2,avededxtrk2Last4Cells);
358  sandbox.SetAvedEdxTrack2Last4Cells(avededxtrk2Last4Cells);
359 
360  fNumuSandFxs->getAveTrackdEdxLast6Cells(track2,avededxtrk2Last6Cells);
361  sandbox.SetAvedEdxTrack2Last6Cells(avededxtrk2Last6Cells);
362 
363  fNumuSandFxs->getAveTrackdEdxLast8Cells(track2,avededxtrk2Last8Cells);
364  sandbox.SetAvedEdxTrack2Last8Cells(avededxtrk2Last8Cells);
365 
366  scattangletrk2 = fNumuSandFxs->getScatt(track2);
367  if(track2->TotalLength() > 0) sandbox.SetScattAngleTrack2(scattangletrk2/track2->TotalLength());
368  //}
369 
370  ///////////////////////////////////
371  ///--------------------------------
372  /// ADD HERE VARS FOR QE PID?
373  ///-------------------------------
374  ///////////////////////////////////
375 
376  // WORKING AREA
377  //art::PtrVector<rb::CellHit> allslchits = slice->AllCells();
378  art::PtrVector<rb::CellHit> alltrk1hits = track1->AllCells();
379  art::PtrVector<rb::CellHit> alltrk2hits = track2->AllCells();
380 
381  fNumuSandFxs->getActivityVtx(alltrk1hits,alltrk2hits,allslchits,actvtxx);
382  slccalE = slice->CalorimetricEnergy();
383 
384  if(slice->CalorimetricEnergy() <= 0){
385  //continue;
386  std::cout<<"Calorimetric energy is zero!!"<<std::endl;
387  }
388  else{
389  sandbox.SetActVtx(actvtxx/slccalE);
390  //std::cout<<"TMVAvars[0]: "<< actvtxx/slccalE <<std::endl;
391  //TMVAvars[0] = actvtx/slccalE;
392  }
393 
394  //----
395  // Another working area
396  // from QEEventFinder/QeFinder_module.cc
397  // get numu energy estimates
398  float trkCCE = -5.0;
399  float RecoMuonE = -5.0;
400  float Reco2ndTrkE = -5.0;
401 
402  if(eAss.isValid()){
403  std::vector<art::Ptr<numue::NumuE> > numuEnergy = eAss.at(i);
404  if(numuEnergy.size() > 0){
405  // get the energy variables
406  trkCCE = numuEnergy[0]->TrkCCE();
407  }//end if valid size
408 
409  art::FindOneP<rb::Energy> trackEnergy(sliceTracks, evt, fEnergyLabel);
410  if(trackEnergy.isValid()){
411  RecoMuonE = trackEnergy.at(bestIdx)->E();
412  }
413  //Reco2ndTrkE = track2->TotalGeV(); // sloppy approach
414  Reco2ndTrkE = track2->CalorimetricEnergy(); // sloppy approach
415 
416  float tmpHadE = trkCCE - RecoMuonE;
417 
418  float offTrkEnFrac = (tmpHadE - Reco2ndTrkE)/trkCCE;
419 
420  sandbox.SetOffTrkFra(offTrkEnFrac);
421 
422  //std::cout<<"TMVAvars[0]: "<< offTrkEnFrac <<std::endl;
423 
424  }// eAss is valid
425  else{
426  std::cout<<"Event Association is not Valid!"<<std::endl;
427  }
428 
429  }//end slice tracks = 2
430 
431  // ----------------------------------
432 
433  // Loop over tracks in slice
434  for(size_t j = 0; j < sliceTracks.size(); ++j){
435  const art::Ptr<rb::Track> track = sliceTracks[j];
436  art::Ptr<remid::ReMId> trackRemid = remidVec.at(j);
437 
438  if(!track->Is3D()) {
439  continue;
440  }
441 
442  idx3d++;
443 
444  if(trackRemid->Value() > bestval) { // best track so far
445  val2 = bestval;
446  Idx2 = bestIdx;
447  bestval = trackRemid->Value();
448  bestIdx = j;
449  }
450  else if(trackRemid->Value() > val2) { // not best, but better second best track
451  val2 = trackRemid->Value();
452  Idx2 = j;
453  }
454 
455  if(trackRemid->Value() > 0.6) nmutrks++;
456 
457  if(bt->HaveTruthInfo()) {
458  art::PtrVector<rb::CellHit> trkhits = track->AllCells();
459  const std::vector<const sim::Particle*> parts = bt->HitsToParticle(trkhits);
460  const std::vector< cheat::TrackIDE > trids = bt->HitsToTrackIDE(trkhits);
461 
462  // check particle that contributes most to track for pi->mu decay
463  if( !parts.empty() && !trids.empty() ) {
464  if(abs(parts[0]->PdgCode())==13 && parts[0]->Process()=="Decay") {
465  double tempE = trids[0].energyFrac * track->TotalGeV();
466  if(tempE > bestpimuE) {
467  bestpimuE = tempE;
468  pimuIdx = idx3d+1;
469  }
470  }
471  }
472 
473  // check particle that contributes second most; track can contain both pi and mu
474  if(parts.size() > 1 && trids.size() > 1) {
475  if(abs(parts[1]->PdgCode())==13 && parts[1]->Process()=="Decay") {
476  double tempE = trids[1].energyFrac * track->TotalGeV();
477  if(tempE > bestpimuE) {
478  bestpimuE = tempE;
479  pimuIdx = idx3d+1;
480  }
481  }
482  }
483  }
484  } //end loop over tracks
485 
486  // ===============================================================
487  //New: Get info on muon like track
488  //if( (bestval > 0.0) && (val2 > 0.0) ){
489  if((bestval > 0.0)){
490 
491  const art::Ptr<rb::Track> mutrack = sliceTracks[bestIdx];
492  art::PtrVector<rb::CellHit> mutrkhits = mutrack->AllCells();
493 
494  //Getting all the hits in the slice that are not on our best track
495  ////rb::Cluster vertexHits = MakeVertexCluster(allslchits,mutrkhits);
496  rb::Cluster vertexCluster;
497  int nhadhits = -5;
498  for (unsigned int sliceHit = 0; sliceHit < allslchits.size(); ++sliceHit){
499 
500  bool match = 0;
501 
502  for (unsigned int trackHit = 0; trackHit < mutrkhits.size(); ++trackHit){
503  match = (allslchits[sliceHit] == mutrkhits[trackHit]);
504  if (match) break;
505  }//End of loop over hits in track
506 
507  if (!match) vertexCluster.Add(allslchits[sliceHit]);
508 
509  }//End of loop over hits in the slice
510 
511  nhadhits = vertexCluster.NCell();
512  //std::cout<<" ==== No. hits in had system: "<<vertexCluster.NCell()<<std::endl;
513 
514  sandbox.SetNHadHits(nhadhits);
515 
516  // New
517  unsigned int minCellsFromEdge = 99999999;
518  //int cellsFromEdge = -1;
519 
520  // TO DO: case when nhadhits is zero?
521  // Loop over vertex cluster
522  for(unsigned int hitIdx = 0; hitIdx < vertexCluster.NCell(); ++hitIdx){
523  const art::Ptr<rb::CellHit>& chit = vertexCluster.Cell(hitIdx);
524  const int planeNum = chit->Plane();
525  const unsigned int cellsFromEdge = std::min((unsigned int)chit->Cell(), geom->Plane(planeNum)->Ncells() - 1 - chit->Cell());
526 
527  if(cellsFromEdge < minCellsFromEdge){
528  minCellsFromEdge = cellsFromEdge;
529  } //end if
530 
531  } //end for
532  //std::cout<<" ==== Had ncells from edge: "<< minCellsFromEdge <<std::endl;
533 
534  sandbox.SetHadCellsFromEdge(minCellsFromEdge);
535 
536  }//end if best val
537  // ===============================================================
538 
539  sandbox.SetNMuTrks(nmutrks);
540 
541  if(bt->HaveTruthInfo()) { // no truth info case should correspond to unfilled
542 
543  if(bestIdx!=-1) {
544 
545  const art::Ptr<rb::Track> track1 = sliceTracks[bestIdx];
546  const std::vector<const sim::Particle*> parts1 = bt->HitsToParticle(track1->AllCells());
547  const std::vector< cheat::TrackIDE > trids1 = bt->HitsToTrackIDE(track1->AllCells());
548 
549  if(!parts1.empty()){
550  sandbox.SetPDGBest(parts1[0]->PdgCode());
551  }
552  }
553 
554  if(Idx2!=-1) {
555 
556  const art::Ptr<rb::Track> track2 = sliceTracks[Idx2];
557  const std::vector<const sim::Particle*> parts2 = bt->HitsToParticle(track2->AllCells());
558  const std::vector< cheat::TrackIDE > trids2 = bt->HitsToTrackIDE(track2->AllCells());
559 
560  if(!parts2.empty()){
561  sandbox.SetPDGSecondBest(parts2[0]->PdgCode());
562  }
563  }
564 
565  int nprotons = 0;
566  if(pimuIdx == -5) pimuIdx = 0; // differentiate unfilled with not found
567  art::PtrVector<rb::CellHit> slchits = slice->AllCells();
568  const std::vector<const sim::Particle*> parts = bt->HitsToParticle(slchits);
569  const std::vector< cheat::TrackIDE > trids = bt->HitsToTrackIDE(slchits);
570  if(!parts.empty() && !trids.empty()) {
571  for(unsigned int i = 0; i < parts.size(); i++) {
572  if(trids.size() > i) {
573  if(parts[i]->PdgCode()==2212 && trids[i].energyFrac * slice->TotalGeV() > 0.05) nprotons++;
574  }
575  }
576  }
577  sandbox.SetNProtons(nprotons);
578  }
579 
580  sandbox.SetPiMuDecay(pimuIdx);
581 
582  outputObjects->push_back(sandbox);
583  util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
584  } //end loop over slices
585 
586  evt.put(std::move(outputObjects));
587  evt.put(std::move(assoc));
588 
589  } //end produce
590 
591 } //end namespace
592 
593 
void SetNHadHits(int nhadhits)
void SetXMaxt(float dxmaxt)
void SetPDGSecondBest(int pdgsecondbest)
back track the reconstruction to the simulation
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.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
std::map< std::string, double > xmax
void beginRun(art::Run &run)
void SetZMaxt(float dzmaxt)
void SetXMint(float dxmint)
void SetAvedEdxTrack1Last4Cells(float avededxtrk1Last4Cells)
unsigned short Plane() const
Definition: CellHit.h:39
void SetYMax2(float dymax2)
Definition: NumuSandObj.cxx:92
void SetAvedEdxTrack2Last4Cells(float avededxtrk2Last4Cells)
void SetOffTrkFra(float offtrk)
proxy for off track energy
geo::View_t View() const
Definition: CellHit.h:41
void getAveTrackdEdx(const art::Ptr< rb::Track > track, float &avededx)
T sqrt(T number)
Definition: d0nt_math.hpp:156
void getAveTrackdEdxLast4Cells(const art::Ptr< rb::Track > track, float &avededxlast4)
Vertical planes which measure X.
Definition: PlaneGeo.h:28
FillSandbox(fhicl::ParameterSet const &pset)
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
float Z() const
Definition: RecoHit.h:38
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
float getScatt(const art::Ptr< rb::Track > track)
void abs(TH1 *hist)
void SetAvedEdxTrack2Last6Cells(float avededxtrk2Last6Cells)
void SetScattAngleTrack1(float scattangletrk1)
void SetAvedEdxTrack2Last8Cells(float avededxtrk2Last8Cells)
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
void SetPDGBest(int pdgbest)
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
double Value() const
Definition: PID.h:22
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void SetXMin2(float dxmin2)
minimum position of hit second farthest from edge in that dimension
Definition: NumuSandObj.cxx:68
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
Definition: Run.h:31
void SetNMuTrks(int nmutrks)
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
void SetPiMuDecay(int pitomu)
Definition: NumuSandObj.cxx:56
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
void SetVtxE20(float vtxE20)
Double_t ymax
Definition: plot.C:25
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
void produce(art::Event &evt)
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
void getActivityVtx(const art::PtrVector< rb::CellHit > track1Hits, const art::PtrVector< rb::CellHit > track2Hits, const art::PtrVector< rb::CellHit > sliceHits, float &actvtx)
void SetYMaxt(float dymaxt)
minimum position of hit closest to edge in that dimension that has more than 100 photoelectons (not p...
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Definition: Cluster.cxx:439
void SetHadCellsFromEdge(int hadcellsfromedge)
void SetYMint(float dymint)
minimum position of hit closest to edge in that dimension that has more than 100 photoelectons (not p...
void SetAvedEdxTrack1Last6Cells(float avededxtrk1Last6Cells)
int evt
void SetAvedEdxTrack1(float avededxtrk1)
float PE() const
Definition: CellHit.h:42
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
Float_t d
Definition: plot.C:236
void getAveTrackdEdxLast8Cells(const art::Ptr< rb::Track > track, float &avededxlast8)
void SetZMint(float dzmint)
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
const double j
Definition: BetheBloch.cxx:29
void SetAvedEdxTrack1Last8Cells(float avededxtrk1Last8Cells)
void SetAvedEdxTrack2(float avededxtrk2)
Vertex location in position and time.
z
Definition: test.py:28
Sandbox to make objects for numu analysis in CAF.
Definition: FillSandboxes.h:7
void SetZMax2(float dzmax2)
Definition: NumuSandObj.cxx:98
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:308
void SetMissingE(float misse)
total E of true particles leaving detector (>1 MeV, not neutron)
Definition: NumuSandObj.cxx:62
OStream cout
Definition: OStream.cxx:6
void SetYMin2(float dymin2)
Definition: NumuSandObj.cxx:74
numusand::NumuSandFxs * fNumuSandFxs
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
std::vector< TrackIDE > HitsToTrackIDE(const std::vector< const rb::CellHit * > &hits, bool useBirksE=false) const
Returns vector of TrackIDE structs contributing to the given collection of hits.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void SetNDropped(int ndropped)
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
float GeV() const
Definition: RecoHit.cxx:69
void getAveTrackdEdxLast6Cells(const art::Ptr< rb::Track > track, float &avededxlast6)
void SetXMax2(float dxmax2)
maximum position of hit second farthest from edge in that dimension
Definition: NumuSandObj.cxx:86
T const * get() const
Definition: Ptr.h:321
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
float X() const
Definition: RecoHit.h:36
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
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
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
void SetActVtx(float startact)
proxy for hadE[all hits not in main track] - calE[2nd reco trk]
Double_t ymin
Definition: plot.C:24
void SetScattAngleTrack2(float scattangletrk2)
bool getMissingE(const art::Ptr< rb::Cluster > slice, float &missE1)
void SetZMin2(float dzmin2)
Definition: NumuSandObj.cxx:80
float Y() const
Definition: RecoHit.h:37
void SetNBad(int nbad)
T min(const caf::Proxy< T > &a, T b)
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
Definition: Prong.h:71
void SetVtxE60(float vtxE60)
void SetNProtons(int nprotons)
void SetVtxE40(float vtxE40)
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.