BremShowerFilter_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file BremShowerFilter.cxx
3 /// \author H. Duyang 12/18/13
4 ////////////////////////////////////////////////////////////////////////
5 // ART includes
14 #include "fhiclcpp/ParameterSet.h"
18 
19 // NOvA includes
20 #include "Geometry/Geometry.h"
21 //#include "Geometry/Geo.h"
23 #include "Calibrator/Calibrator.h"
24 #include "MCCheater/BackTracker.h"
27 #include "RecoBase/Track.h"
28 #include "RecoBase/Cluster.h"
29 #include "RecoBase/CellHit.h"
30 #include "RecoBase/RecoHit.h"
31 #include "RawData/DAQHeader.h"
32 #include "RawData/RawDigit.h"
33 #include "Utilities/func/MathUtil.h"
34 #include "NovaDAQConventions/DAQConventions.h"
35 #include "DAQDataFormats/NanoSliceVersionConvention.h"
36 
37 //root includes
38 #include <vector>
39 #include "TH1F.h"
40 #include "TH2F.h"
41 #include "TH3F.h"
42 #include "TTree.h"
43 #include "TLorentzVector.h"
44 
45 namespace bsf {
46 
48 
49  public:
50  explicit BremShowerFilter(fhicl::ParameterSet const& pset);
51  virtual ~BremShowerFilter();
52 
53  virtual bool filter(art::Event & evt);
54  void reconfigure(const fhicl::ParameterSet& pset);
55  void beginJob();
57  void findShowerByTruth(art::Ptr<rb::Track> &trk, int &startPlane, int &endPlane);
58  bool inFiducial(double x, double y, double z);
61 
62  private:
63  bool dokeep;
64  //parameters here
65 
69 
70  std::string fTrackLabel; ///< label of module making fls hits
71 
74 
75  //find shower by truth or by reco
76  bool fShowerByTruth;//if true, look for only DIF using truth
77  bool fShowerByReco;//if true look for all shower using reco
78  bool fFindDiF;//if true, look for only DIF using reco
79  float fAngleCut;//cosZ cut on muon track
80  float fPlaneCut;//how many plane out of all planes we require?
81  float fEnergyCut;//requre a minimum shower energy
82  float fFiducialX;
83  float fFiducialY;
84  float fFiducialZ;
85  int fPlaneMin;
86  int fPlaneMax;
87 
88  //make a tree to save variables for analysis
89  TTree *fTree;
90 
91  int run;
92  int subrun;
93  int event;
94  int particleId;//this is the true particle id
95  int trackId;//this is the reco track id
96  //int process;
97  //int decayPlane;
99  int endPlane;
100  //int nplane;
101  double cosz;
102  double eshower;
103  double planeDedx[2];
104  double planeDedxMR[2];
105  double cellDedx[2];
106  double cellDedxMR[2];
108 
109 
110  //double E_mu_decay;
111  //double ADCTrack[1000][500];
112  //double edepTrack[1000][500];
113  //double ADCSlice[1000][500];
114  //double edepSlice[1000][500];
115  //double ADCShower[1000][500];
116  //double edepShower[1000][500];
117  //double edepEl[1000][500];
118  //double edepMu[1000][500];
119 
120  //TLorentzVector *v4; // particle vertex
121  //TLorentzVector *p4Mu; // particle 4 momentum
122  //TLorentzVector *p4El; // mother 4 momentum
123  //TVector3 *dir;
124  TVector3 *vtx= new TVector3(0,0,0);
125 
126  //hists here
127  //TH1F *fHEmu[5];
128  //TH1F *fHEndPlane;
129  //TH1F *fHEndCell;
130  TH1F *fHNTrk;
131  TH1F *fHCosZ;
132  TH1F *fHCosZCut[5];
135  TH1F *fHNPlane;
136  TH1F *fHEShower;
137  //TH3F *fHTrackStart;
138  //TH3F *fHTrackEnd;
140  //TH2F *fHEndPoint;
141  TH1F *fHEdep;
142  TH1F *fHADC;
143  TH1F *fHPE;
144  TH1F *fHPECorr;
145 
146  TH2F *fHHitMap;
147  TH2F *fHHitMapRW;
148 
149  };
150 }
151 
152 namespace bsf {
153  //......................................................................
154 
156  : fTrackToken(consumes<std::vector<rb::Track>>(pset.get<std::string>("TrackLabel"))),
157  fClusterToken(consumes<std::vector<rb::Cluster>>(pset.get<std::string>("ClusterLabel"))),
158  fCellHitToken(consumes<std::vector<rb::CellHit>>(pset.get<std::string>("CellHitLabel")))
159  //v4(), p4Mu(), p4El()
160  {
161  reconfigure(pset);
162  //produces<std::vector<dm::BremShowerInfo> >();
163  produces< std::vector<rb::Cluster> >();
164  produces<std::vector<rawdata::RawDigit> >();
165 
166  //produces< std::vector<rb::Track> >();
167  }
168 
169  //......................................................................
170 
172 
173  //......................................................................
174 
175 
178 
179  fTree = tfs->make<TTree>("brem", "brem");
180  fTree->Branch("run", &run, "run/I");
181  fTree->Branch("subrun", &subrun, "subrun/I");
182  fTree->Branch("event", &event, "event/I");
183  fTree->Branch("particleId", &particleId, "particleId/I");
184  fTree->Branch("trackId", &trackId, "trackId/I");
185  //fTree->Branch("decayPlane", &decayPlane, "decayPlane/I");
186  //fTree->Branch("process", &process, "process/I");
187  fTree->Branch("startPlane", &startPlane, "startPlane/I");
188  fTree->Branch("endPlane", &endPlane, "endPlane/I");
189  //fTree->Branch("nPlane", &nPlane, "nPlane/I");
190  fTree->Branch("cosz", &cosz, "cosz/D");
191  fTree->Branch("eshower", &eshower, "eshower/D");
192  fTree->Branch("vtx", "TVector3", &vtx);
193  fTree->Branch("planeDedx", &planeDedx, "planeDedx[2]/D");
194  fTree->Branch("planeDedxMR", &planeDedxMR, "planeDedxMR[2]/D");
195  fTree->Branch("cellDedx", &cellDedx, "cellDedx[2]/D");
196  fTree->Branch("cellDedxMR", &cellDedxMR, "cellDedxMR[2]/D");
197  //fTree->Branch("E_mu_decay", &E_mu_decay, "E_mu_decay/D");
198  //fTree->Branch("v4", "TLorentzVector", &v4);
199  //fTree->Branch("p4El", "TLorentzVector", &p4El);
200  //fTree->Branch("p4Mu", "TLorentzVector", &p4Mu);
201  //fTree->Branch("dir", "TVector3", &dir);
202 
203  //fTree->Branch("edepMu", &edepMu, "edepMu[1000][500]/D");
204  //fTree->Branch("edepEl", &edepEl, "edepEl[1000][500]/D");
205  //fTree->Branch("edepTrack", &edepTrack, "edepTrack[1000][500]/D");
206  //fTrackTree->Branch("ADCTrack", &ADCTrack, "ADCTrack[1000][500]/D");
207  //fTree->Branch("edepSlice", &edepSlice, "edepSlice[1000][500]/D");
208  //fTrackTree->Branch("ADCSlice", &ADCSlice, "ADCSlice[1000][500]/D");
209  //fTree->Branch("edepShower", &edepShower, "edepShower[1000][500]/D");
210  //fTrackTree->Branch("ADCShower", &ADCShower, "ADCShower[1000][500]/D");
211 
212  fHNTrk = tfs->make<TH1F>("hNTrk", "", 100, 0, 100);
213 
214  fHCosZ = tfs->make<TH1F>("hCosZ", "", 200, -1, 1);
215  for(int i=0; i<5; i++)
216  fHCosZCut[i] = tfs->make<TH1F>(Form("hCosZCut%d", i), "", 200, -1, 1);
217  fHNPlane = tfs->make<TH1F>("hNPlane", "", 1000, 0, 1000);
218  fHStartContainment = tfs->make<TH1F>("hStartContainment", "", 2, 0, 2);
219  fHEndContainment = tfs->make<TH1F>("hEndContainment", "", 2, 0, 2);
220  fHEShower = tfs->make<TH1F>("hEShower", "", 1000, 0, 10);
221  fHEdep = tfs->make<TH1F>("hEdep", "", 1000, 0, 1);
222  fHADC = tfs->make<TH1F>("hADC", "", 1000, 0, 1000);
223  fHPE = tfs->make<TH1F>("hPE", "", 1000, 0, 1000);
224  fHPECorr = tfs->make<TH1F>("hPECorr", "", 1000, 0, 1000);
225  //fHTrackStart = tfs->make<TH3F>("hTrackStart", "", 200, -1000, 1000, 200, -1000, 1000, 6000, 0, 6000);
226  //fHTrackEnd = tfs->make<TH3F>("hTrackEnd", "", 200, -1000, 100, 200, -1000, 1000, 6000, 0, 6000);
227  fHTrackStartZ = tfs->make<TH1F>("hTrackStartZ", "", 6000, 0, 6000);
228  fHHitMap = tfs->make<TH2F>("hHitMap", "", 200, 0, 200, 400, 0, 400);
229  fHHitMapRW = tfs->make<TH2F>("hHitMapRW", "", 200, 0, 200, 400, 0, 400);
230 
231  }
232 
233  //......................................................................
234 
236  {
237  fTrackLabel = pset.get< std::string >("TrackLabel");
238  fWriteParentSliceMode = pset.get< bool >("WriteParentSliceMode");
239  fApplyFilter = pset.get< bool >("ApplyFilter");
240  fShowerByTruth = pset.get< bool >("ShowerByTruth", "false");
241  fShowerByReco = pset.get< bool >("ShowerByReco", "false");
242  fFindDiF = pset.get< bool >("FindDiF", "false");
243  fAngleCut = pset.get< float >("AngleCut");
244  fPlaneCut = pset.get< float >("PlaneCut");
245  fEnergyCut = pset.get< float >("EnergyCut");
246  fFiducialX = pset.get< float >("FiducialX");
247  fFiducialY = pset.get< float >("FiducialY");
248  fFiducialZ = pset.get< float >("FiducialZ");
249  fPlaneMin = pset.get< int >("PlaneMin");
250  fPlaneMax = pset.get< int >("PlaneMax");
251  fDCSDistance = pset.get< int >("DCSDistance");
252  }
253 
254  //......................................................................
256  {
259 
260  //Look at the truth
261  //Get the particle navigator
262  //art::ServiceHandle<cheat::BackTracker> bt;
263  //const sim::ParticleNavigator &pnav = bt->ParticleNavigator();
264 
265  //use reco base
267  evt.getByToken(fTrackToken, trackcol);
268 
270  evt.getByToken(fCellHitToken, cellhit);
271 
272  // Get the slices in the event
274  evt.getByToken(fClusterToken, slicecol);
275 
276  art::FindManyP<rb::Track> fmt(slicecol, evt, fTrackLabel);
277 
278  std::vector< rawdata::RawDigit> newDigits;
279  newDigits.reserve(cellhit->size());
280  for(const rawdata::RawDigit& digit: *cellhit) newDigits.push_back(digit);
281 
282  //std::unique_ptr<std::vector<dm::BremShowerInfo> >
283  //BremShower(new std::vector<dm::BremShowerInfo>);
284 
285  //and the rawdigit for the showers
286  std::unique_ptr< std::vector<rawdata::RawDigit> > showerDigits(new std::vector<rawdata::RawDigit>);
287 
288  //a new cluster to hold the showers
289  std::unique_ptr< std::vector<rb::Cluster> > showercol(new std::vector<rb::Cluster>);
290  //std::unique_ptr< std::vector<rb::Track> > etrackcol(new std::vector<rb::Track>);
291 
292  dokeep = false; //indicate to save the event or not
293 
294  run = evt.run();
295  subrun = evt.subRun();
296  event = evt.id().event();
297 
298  fHNTrk->Fill(trackcol->size());
299 
300  //prepare a vector of hits that needs to be removed
301  std::vector<rb::WeightedHit > trackHits;//while I call this track hits, it is actually all hits except those belong to ther shower.
302 
303  //loop over slices here
304  for(unsigned int iSlice = 0; iSlice < slicecol->size(); ++iSlice) {
305 
306  art::Ptr<rb::Cluster> slc(slicecol,iSlice);
307 
308  bool goodSlice = true;
309  if (slc->NCell() <= 0) continue;
310  if(slc->MeanTNS () <25000 || slc->MeanTNS () > 475000) {goodSlice = false;} // rejecting the trigger window on edges
311 
312  std::vector< art::Ptr<rb::Track> > tracks = fmt.at(iSlice);
313 
314  //is there a shower on this track?
315  //if yes, where it starts and where it ends?
316  //assume there is no more than one shower per track for now
317 
318  startPlane = -1;
319  endPlane = -1;
320  trackId = iSlice; //what saving here is actually the slice where the shower from
321 
322  bool isShower = false;
323  bool isCandTrack = false;
324  //bool isBrem = false;
325  bool isDiF = false;
326 
327  if(!slc->IsNoise() && tracks.size()==1 && tracks[0]->Is3D() && goodSlice){//if there is no track, then no shower at all
328 
329  bool startContained = false;
330  bool endContained = false;
331 
332  //first decide is this track an candidate track?
333  if(int(tracks[0]->MaxPlane()) > fPlaneMin && int(tracks[0]->MinPlane()) < fPlaneMax){//well this means the track passes the active area of the detector
334  startContained = trackStartContained(tracks[0]);
335  endContained = trackEndContained(tracks[0]);
336 
337  fHStartContainment->Fill(startContained);
338  fHEndContainment->Fill(endContained);
339 
340  //fHTrackStart->Fill(tracks[0]->Start().X(), tracks[0]->Start().Y(), tracks[0]->Start().Z());
341  fHTrackStartZ->Fill(tracks[0]->Start().Z());
342 
343  int nPlaneTrack = 0;
344 
345  if(int(tracks[0]->MaxPlane())>fPlaneMax) nPlaneTrack = fPlaneMax - int(tracks[0]->MinPlane());
346  else if(int(tracks[0]->MinPlane())<fPlaneMin) nPlaneTrack = int(tracks[0]->MaxPlane()) - fPlaneMin;
347  else nPlaneTrack = tracks[0]->MaxPlane() - tracks[0]->MinPlane();
348 
349  cosz = tracks[0]->Dir().Z();
350 
351  if(startContained) isCandTrack = false;//#1 cut on containment
352  else{
353  fHCosZ->Fill(tracks[0]->Dir().Z());
354  //if(fabs(tracks[0]->Dir().Z())<fAngleCut) isCandTrack = false;//#2 cut on angle
355  if(tracks[0]->Dir().Z()<fAngleCut) isCandTrack = false;//#2 cut on angle, keep only downstream tracks
356  else {
357  fHNPlane->Fill(nPlaneTrack); //#3 cut on nPlane
358  if(nPlaneTrack<fPlaneCut) isCandTrack = false;
359  else{
360  isCandTrack = true;
361  fHCosZCut[0]->Fill(tracks[0]->Dir().Z());
362  }
363  }
364  }
365 
366  }//end if track passes the activ area of the detector
367 
368  if(isCandTrack){//we found the candidate track!
369  if(fShowerByReco){//means to find both brem and DiF
370 
371  findShowerByReco(tracks[0], startPlane, endPlane);
372  if(startPlane>0 && endPlane>0){
373  fHCosZCut[1]->Fill(tracks[0]->Dir().Z());
374  if(startPlane>fPlaneMin && endPlane<fPlaneMax){//select shower contained in the active area of the detector
375  fHCosZCut[2]->Fill(tracks[0]->Dir().Z());
376  isShower = true;
377  //if(endContained)
378  //isDiF = true;
379  }
380 
381  }
382  }//end if fShowerByReco
383 
384  //if(fFindDiF){//means to find only DiF
385  //findShowerByReco(tracks[0], startPlane, endPlane);
386  //if(!startContained && endContained && startPlane>0 && endPlane<0) {isShower = true; isDiF = true;}
387  //}
388 
389  if(fShowerByTruth){
390  findShowerByTruth(tracks[0], startPlane, endPlane);
391  if(!startContained && endContained && startPlane>0 && endPlane>0) {isShower = true; isDiF = true;}
392  }
393  }//end if is candidate track
394 
395  }//end if track size is 1
396 
397  //the following part fills a vector of weighted hits trackHits to be removed
398  //consider to use a function for this
399  //getTrackHits(tracks[0], startPlane, endPlane, trackHits)
400 
401  if(isShower){//we found the shower!
402  //in this case, keep the hits in shower region, remove the rest in the slice
403  //also save the shower to a new cluster object
404  //bool muonRemove = false;
405  //if(endPlane>0) muonRemove=true;//we call it brem. This criteria is not a good one of course
406 
407  //prepare an cluster object to save the shower hits
408  rb::Cluster shower;
409 
410  //shower region should be taken care of first because we want to remove some hit while keep others
411  //we need a loop over planes in shower region!
412 
413  int minPlane = startPlane;
414  int maxPlane;
415  if(endPlane>startPlane) maxPlane = endPlane;
416  else maxPlane = tracks[0]->MaxPlane();
417 
418 
419  TVector3 startPoint;
420  geom->Plane( startPlane )->Cell( 10 )->GetCenter(startPoint);
421 
422  //try to define a shower vertex here
423  //loop over trajectory points, find the one on the startPlane, call it the vertex
424  //note that it is only true for a down stream shower.
425  int nTrajPt = tracks[0]->NTrajectoryPoints();
426 
427  for( int iTraj = 0 ; iTraj < nTrajPt ; iTraj++ ){
428  //this is to find the vertex
429  if( startPoint.Z()-3.3 < tracks[0]->TrajectoryPoint(iTraj).Z() &&
430  startPoint.Z()+3.3 > tracks[0]->TrajectoryPoint(iTraj).Z()){
431  vtx->SetXYZ(tracks[0]->TrajectoryPoint(iTraj).X(),
432  tracks[0]->TrajectoryPoint(iTraj).Y(),
433  tracks[0]->TrajectoryPoint(iTraj).Z());
434 
435  }
436 
437 
438  }//end loop over trajectory points
439 
440  //loop over planes in shower region
441  for(int iPlane=minPlane; iPlane<=maxPlane; iPlane++){
442  std::vector< art::Ptr<rb::CellHit> > planeHits;//prepare a vector of cell hits on this plane
443 
444  //the cell on this plane closest, and the 2nd closest to track,
445  unsigned int cellClose1=0;
446  unsigned int cellClose2=0;
447  double cellCloseDist1=10000000;
448  double cellCloseDist2=10000000;
449  double edepClose1 = 0;
450  double edepClose2 = 0;
451 
452  //int trajpoint = -1;
453 
454  //loop over trajectory points
455  //for( int iTraj = 1; iTraj < nTrajPt ; iTraj++){
456 
457  //const geo::CellUniqueId cid = geom->CellId(tracks[0]->TrajectoryPoint(iTraj).X(),
458  // tracks[0]->TrajectoryPoint(iTraj).Y(),
459  // tracks[0]->TrajectoryPoint(iTraj).Z(),
460  // 1,1,0, 2);
461  //int trajPlane = -1;
462 
463  //geom->IdToPlane( cid, &trajPlane);
464 
465  //if(trajPlane==iPlane){
466  //trajpoint = iTraj;
467  //break;
468  //}
469  //}//end loop over trajectory points
470 
471  //double trajXY = 0;
472  //if(trajpoint>0){
473  //if(geom->Plane(iPlane)->View() == geo::kX)
474  // trajXY = tracks[0]->TrajectoryPoint(trajpoint).X();
475  //else
476  // trajXY = tracks[0]->TrajectoryPoint(trajpoint).Y();
477  //}
478  //else continue;//not trajpoint found in this plane; use the second method. (I don't think this will happen)
479 
480  //this is the energy a muon *should* deposit on this plane
481  double mipDeDx = 0.00157;
482  //double mipDeDx = 0.00135;
483  double mipPlane = 2*mipDeDx*geom->Plane(iPlane)->Cell(0)->HalfD()/fabs(tracks[0]->Dir().Z());
484 
485  //then loop over track hits, for those in shower region, fill a vector with hits on this plane
486  //
487  for(unsigned int iCell=0; iCell<tracks[0]->NCell(); ++iCell){
488  const rb::RecoHit rhit = tracks[0]->RecoHit(tracks[0]->Cell(iCell));
489  const rb::CellHit *chit = tracks[0]->Cell(iCell).get();
490 
491  if(!rhit.IsCalibrated()) continue;
492  if(chit->Plane() != iPlane) continue;
493  planeHits.push_back(tracks[0]->Cell(iCell));
494 
495  //decide what is the distance of the hits to track, if too far away, just remove it(will do)
496  double xyz[3];
497  geom->Plane( chit->Plane() )->Cell( chit->Cell() )->GetCenter(xyz);
498 
499  //well this is a stupid way of calculating the distance
500  double cellXY, trackStartXY;
501  double cellZ = xyz[2];
502  double trackStartZ = tracks[0]->Start().Z();
503  double slope = tracks[0]->Dir().X()/tracks[0]->Dir().Z();
504 
505  if(geom->Plane(chit->Plane())->View() == geo::kX){
506  cellXY = xyz[0];
507  trackStartXY = tracks[0]->Start().X();
508  slope = tracks[0]->Dir().X()/tracks[0]->Dir().Z();
509  }
510  else{
511  cellXY = xyz[1];
512  trackStartXY = tracks[0]->Start().Y();
513  slope = tracks[0]->Dir().Y()/tracks[0]->Dir().Z();
514  }
515 
516  double cellDist = fabs(cellXY - trackStartXY - slope*(cellZ - trackStartZ));
517  //double cellDist1 = fabs(cellXY - trajXY);
518  //sort plane hits according to the dist form track.
519 
520  if(cellDist<cellCloseDist1){
521  cellCloseDist2 = cellCloseDist1;
522  cellClose2 = cellClose1;
523  edepClose2 = edepClose1;
524  cellCloseDist1 = cellDist;
525  cellClose1 = chit->Cell();
526  edepClose1 = rhit.GeV();
527  }
528  else if(cellDist<cellCloseDist2){
529  cellCloseDist2 = cellDist;
530  cellClose2 = chit->Cell();
531  edepClose2 = rhit.GeV();
532  }
533  }//end loop over track hits
534 
535  //loop over plane hits... try to define a weight and fill weighted hits
536  for(unsigned int iHit = 0; iHit < planeHits.size(); ++iHit){
537 
538  art::Ptr<rb::CellHit> chit = planeHits[iHit];
539  //rb::RecoHit rhit = tracks[0]->RecoHit(chit->Cell());
540 
541  double weight = 0;//the default should be zero, mean don't remove
542  if(!isDiF){//only for brem we need to remove the muon in shower region
543 
544  if(chit->Cell()==cellClose1){//the cell closest to track
545  if(edepClose1>mipPlane){//if this cell alone has energy > mip
546  //reset the energy to be a mip
547  weight = mipPlane/edepClose1;
548  if(weight<0 || weight>1) {
549  mf::LogWarning("BremShowerFilter") <<"weight = "<<weight<<" "<<mipPlane<<" "<<edepClose1<<" "<<edepClose2;
550  }
551  }
552  else //if the cell has energy <= mip
553  weight = 1;//fill this one with weight=1, and remove it later
554  }
555  else if(chit->Cell()==cellClose2){//the 2nd closest cell
556  //if the sum of the first and the second hits larger than a mip, reset the energy of the second to be what left
557  if(edepClose1<mipPlane){
558  if(edepClose1+edepClose2>mipPlane){
559  weight = (-edepClose1 + mipPlane)/edepClose2;
560  if(weight<0 || weight>1) {
561  mf::LogWarning("BremShowerFilter") <<"weight = "<<weight<<" "<<mipPlane<<" "<<edepClose1<<" "<<edepClose2;
562  }
563  }
564  else //if the two closest cells still don't have enough energy to make the mip, remove them all. leave the rest.
565  weight = 1;
566  //otherwise do nothing, weight = 1
567  }
568  // else continue;//if the closest aleady has > mip energy, don't touch the second
569  }
570  else weight=0;//keep all other cell hits, means weight should be 0!
571  //else continue;
572  }//endl if !isDiF
573  else weight = 0;//if it is dif then don't remove anything
574 
575  if(weight<0 || weight>1) {
576  mf::LogWarning("BremShowerFilter") <<"weight = "<<weight;
577  }
578 
579  if(weight>0){//save trackhits
580  rb::WeightedHit newHit(chit, weight);
581  trackHits.push_back(newHit);
582  }
583  else if((1-weight)>0)//save shower hits
584  shower.Add(chit, (1-weight));
585 
586  }//endl loop over plane hits
587  }//end loop over planes
588 
589  //loop over track hits again, this time for removing those outside shower region (I know this method sounds stupid)
590  //4/11/2014: change to slice hits
591  for(unsigned int iCell=0; iCell<slc->NCell(); ++iCell){
592  //art::Ptr<rb::CellHit> chit = slices[i]->Cell(iCell);
593  const rb::CellHit *chit = slc->Cell(iCell).get();
594  int plane = chit->Plane();
595 
596  if(plane<minPlane || plane>maxPlane){//not in shower region, simplay remve the hits
597  rb::WeightedHit newHit(slc->Cell(iCell), 1);
598  trackHits.push_back(newHit);
599  }
600 
601  }//end loop other track hits
602 
603  if(!fWriteParentSliceMode) showercol->push_back(shower);
604 
605  shower.Clear();
606  dokeep = true;//keep this event as long as there is shower
607  fTree->Fill();
608 
609  if(fWriteParentSliceMode) showercol->push_back(*slc);
610  }//end if shower found
611  else{//no shower found. remove all hits in this slice
612  //loop over slice hits. For now we just use the slice with the same index
613  //if(iTrack<slicecol->size())
614  for(unsigned int iCell=0; iCell<slc->NCell(); ++iCell){
615  //art::Ptr<rb::CellHit> chit = slices[i]->Cell(iCell);
616  //const rb::RecoHit rhit = slices[i]->RecoHit(slices[j]->Cell(iCell));
617  rb::WeightedHit newHit(slc->Cell(iCell), 1);
618  trackHits.push_back(newHit);
619  }//endl loop over slice hits
620 
621  }//end else (no shower found)
622 
623  }//end loop over slices
624 
625  //maybe it is faster to remove the trackhits all together here, instead in each track loop
626  //-----------------------------------------------------------------------
627  //prepare all digits in this event
628  //now we want to clean up the track hits
629  //if(!RemoveMuon)
630 
631  if(dokeep){//only when the event has shower we want to remove the non-shower digits. otherwise the event will not ke saved anyway
632  std::map<std::tuple<int, int, int>, rb::WeightedHit> trkMap;
633  for(int iHit = 0; iHit<(int)trackHits.size(); iHit++ ){
634  const rb::WeightedHit& h = trackHits[iHit];
635  trkMap.emplace(std::make_tuple(h.hit->Plane(), h.hit->Cell(), h.hit->TDC()), h);
636  }
637 
638  for(int iDigit = 0; iDigit<(int)newDigits.size(); iDigit++ ){
639 
640  const rawdata::RawDigit& rd = newDigits[iDigit];
641 
642  const rb::CellHit cHit = cal->MakeCellHit( &rd );
643 
644  auto it = trkMap.find(std::make_tuple(cHit.Plane(), cHit.Cell(), cHit.TDC()));
645  if(it == trkMap.end()){
646  showerDigits->push_back(newDigits[iDigit]);
647  continue;
648  }
649 
650  const float muonWeight = it->second.weight;
651 
652  if( muonWeight == 1){
653  // Nothing, hit is not used in output
654  }
655  else{
656  //Mulitpoint readout begins. Copied from MRCC
657  // newDigits[iDigit].SetADC(0, (1-muonWeight)*cHit.ADC(0) );
658 
659  if( newDigits[iDigit].IsMC() && newDigits[iDigit].Version() == 0 && newDigits[iDigit].NADC() == 3){
660  newDigits[iDigit].SetADC(0, (1-muonWeight)*cHit.ADC(0) );
661  newDigits[iDigit].SetADC(1, (1-muonWeight)*cHit.ADC(1));
662  newDigits[iDigit].SetADC(2, (1-muonWeight)*cHit.ADC(2));
663  }
664  else if(newDigits[iDigit].Version() == 0 )
665  newDigits[iDigit].SetADC(0, (1-muonWeight)*cHit.ADC(0) );
666  else{
667 
669  unsigned int nPretrig = conv.getNPretriggeredSamples(newDigits[iDigit].Version());
670  int16_t baseline = newDigits[iDigit].ADC(nPretrig - fDCSDistance );
671  int nadc = newDigits[iDigit].NADC();
672 
673  for( int iadc = 0; iadc < nadc; ++iadc){
674  int16_t newadc = ((newDigits[iDigit].ADC(iadc) - baseline )*( 1-muonWeight) ) + baseline;
675  newDigits[iDigit].SetADC( iadc, newadc);
676  }
677  }
678 
679  showerDigits->push_back(newDigits[iDigit]);
680  }
681  //Mulitpoint readout ends.
682  // newDigits[iDigit].SetADC(0, (1-muonWeight)*cHit.ADC(0) );
683  // newDigits[iDigit].SetADC(1, (1-muonWeight)*cHit.ADC(1));
684  // newDigits[iDigit].SetADC(2, (1-muonWeight)*cHit.ADC(2));
685  }// End loop over digits.
686 
687 
688  //evt.put(std::move(etrackcol));
689  }//end if dokeep
690 
691  evt.put(std::move(showerDigits));
692  evt.put(std::move(showercol));
693 
694  if(fApplyFilter) return dokeep; else return true;
695 
696  }//end filter
697 
698 }
699 
700 bool bsf::BremShowerFilter::inFiducial(double x, double y, double z){
702 
703  bool isContained = false;
704 
705  double xyzMax[3], xyzMin[3];
706 
707  //get the position of the boundary plane
708  geom->Plane(fPlaneMax)->Cell(0)->GetCenter(xyzMax);
709  geom->Plane(fPlaneMin)->Cell(0)->GetCenter(xyzMin);
710 
711  //now get the boundary
712  double zMax = xyzMax[2] + geom->Plane(fPlaneMax)->Cell(0)->HalfD();
713  double zMin = xyzMin[2] - geom->Plane(fPlaneMin)->Cell(0)->HalfD();
714 
715  if(x < (geom->DetHalfWidth() - fFiducialX) && x > (-geom->DetHalfWidth() + fFiducialX) &&
716  y < (geom->DetHalfHeight() - fFiducialY) && y > (-geom->DetHalfHeight() + fFiducialY) &&
717  //z < (geom->DetLength() - fFiducialZ) && z > fFiducialZ
718  z < (zMax - fFiducialZ) && z > zMin + fFiducialZ
719  )
720 
721  isContained = true;
722 
723  return isContained;
724 
725 }
726 
728 
729  bool isContained = inFiducial(trk->Start().X(), trk->Start().Y(), trk->Start().Z());
730 
731  return isContained;
732 }
733 
735  //art::ServiceHandle<geo::Geometry> geom;
736  //loop over track hits, prepare eplane for future use
737  double ymin = 10000;
738 
739  bool isContained = false;
740 
741  /*
742 
743  for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){
744 
745  const rb::CellHit *chit = trk->Cell(iCell).get();
746 
747  double xyz[3] = {0.};
748  geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
749 
750  if(xyz[2]<ymin){
751  isContained = geom->isInsideFiducialVolume(xyz[0], xyz[1], xyz[2]);
752  ymin = xyz[2];
753  }
754  }//end loop over track hits
755  */
756 
757  int nTrajPt = trk->NTrajectoryPoints();
758 
759  for(int iTrajPt=0; iTrajPt<nTrajPt; iTrajPt++){
760  double xTrajpt = trk->TrajectoryPoint(iTrajPt).X();
761  double yTrajpt = trk->TrajectoryPoint(iTrajPt).Y();
762  double zTrajpt = trk->TrajectoryPoint(iTrajPt).Z();
763 
764  if(yTrajpt<ymin){
765  isContained = inFiducial(xTrajpt, yTrajpt, zTrajpt);
766  ymin = yTrajpt;
767  }
768  }
769 
770  return isContained;
771 }
772 
773 //find shower by reco info. the input is rb::Track, output is startPlane, endPlane.
774 //if no shower, startPlane=endPlane=-1
776  //for now, simply cut off tracks too short, or to steep.
777 
779 
780  int minPlane = trk->MinPlane();
781  int maxPlane = trk->MaxPlane();
782  //int nPlane = geom->NPlanes();
783  //fHCosZ->Fill(trk->Dir().Z());
784  //if(fabs(trk->Dir().Z())<fAngleCut) return;
785  //fHNPlane->Fill(fabs(maxPlane-minPlane));
786  //if(fabs(maxPlane-minPlane)<fPlaneCut) return;
787 
788 
789  if(minPlane<fPlaneMin) minPlane = fPlaneMin;
790  if(maxPlane>fPlaneMax) maxPlane = fPlaneMax;
791 
792  std::map<unsigned int, double> ePlane;
793  std::map<unsigned int, double> ePlaneAve;
794 
795  ePlane.clear();
796 
797  //loop over track hits, prepare eplane for future use
798  for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){
799  const rb::CellHit* cHit = trk->Cell(iCell).get();
800  const rb::RecoHit rHit = trk->RecoHit(trk->Cell(iCell));
801 
802  if(cHit->Plane()<fPlaneMin || cHit->Plane()>fPlaneMax) continue;
803 
804  if( rHit.IsCalibrated() ){
805  //ADCTrack[chit->Plane()][chit->Cell()] += chit->ADC(0);
806  //edepTrack[chit->Plane()][chit->Cell()] += rhit.GeV();
807 
808  ePlane[cHit->Plane()] += rHit.GeV();
809  fHEdep->Fill(rHit.GeV());
810  fHADC->Fill(cHit->ADC(0));
811  fHPE->Fill(cHit->PE());
812  fHPECorr->Fill(rHit.PECorr());
813  }
814  }//end loop over cells
815 
816  //start a loop over track planes
817  for(int iPlane=minPlane; iPlane<=maxPlane; iPlane++){
818 
819  if(iPlane == minPlane) ePlaneAve[iPlane] = 0.5*(ePlane[iPlane] + ePlane[iPlane+1]);
820  else if(iPlane == maxPlane) ePlaneAve[iPlane] = 0.5*(ePlane[iPlane-1] + ePlane[iPlane]);
821  else ePlaneAve[iPlane] = (ePlane[iPlane-1] + ePlane[iPlane] + ePlane[iPlane+1])/3;//the average of 3 planes
822  }//end loop over track planes
823 
824  bool isShower = false;
825  bool showerStart = false;
826  bool showerEnd = false;
827 
828  const int nPlaneStart = 5;
829  const int nPlaneEnd = 5;
830 
831  eshower = 0;
832 
833  //start a loop over track planes
834  for(int iPlane=minPlane; iPlane<=maxPlane; iPlane++){
835  double mipDeDx = 0.00157;
836  double mipPlane = 2*mipDeDx*geom->Plane(iPlane)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
837 
838  if(!isShower){//if it's not already in a shower, try to find one
839  //5 planes in a row with energy > 2*mip, call it shower start
840  if(iPlane<=(maxPlane-nPlaneStart))
841  for(int k=0; k<nPlaneStart; ++k){
842  if(k==0) showerStart= (ePlane[iPlane+k]>2*mipPlane);//use eplane or eplaneAve?
843  else showerStart *= (ePlane[iPlane+k]>2*mipPlane);
844  if(!showerStart) break;
845  }
846  else showerStart = false;
847 
848  if(showerStart){//found the shower!
849  isShower = true;//
850  startPlane = iPlane;
851  }//end if shower start
852  }//end if not shower
853  else{//if in the shower region find where it ends
854  for(int k=0; k<nPlaneEnd; k++)
855  if((iPlane+k)<=maxPlane){
856 
857  if(k==0) showerEnd = (ePlaneAve[iPlane+k]<1.5*mipPlane && ePlaneAve[iPlane+k]>0.5*mipPlane);
858  else showerEnd *= (ePlaneAve[iPlane+k]<1.5*mipPlane && ePlaneAve[iPlane+k]>0.5*mipPlane);
859  if(!showerEnd) break;
860  }
861  if(showerEnd){
862  isShower = false;
863  endPlane = iPlane;
864  }
865 
866  }//end else
867 
868  }//end loop over track planes
869 
870  //ok here is another loop over track planes...
871  if(startPlane>0 && endPlane>0 && startPlane>fPlaneMin && endPlane<fPlaneMax){
872  for(int iPlane = startPlane; iPlane <= endPlane; iPlane++) {
873  double mipDeDx = 0.00157;
874  double mipPlane = 2*mipDeDx*geom->Plane(iPlane)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
875  eshower += ePlane[iPlane] - mipPlane;
876  }
877 
878  //nPlane = endPlane - startPlane;
879 
880  fHEShower->Fill(eshower);
881 
882  //if(eshower>0.5) fTree->Fill();
883  }
884 
885 
886  if(eshower < fEnergyCut) startPlane = -1;
887 
888  /*int startPlane;
889  if(trk->Dir().Z()>0) startPlane = minPlane;//going right
890  else startPlane = maxPlane;//going left
891 
892  int iPlane = startPlane;
893  while(iPlane!=endPlane){
894 
895 
896  if(trk->Dir().Z()>0) iPlane++;
897  else iPlane--;
898  }
899  */
900  return;
901 }//end findShowerByReco
902 
903 /*
904 std::vector<rb::WeightedHit> getTrackHits(const rb::Track &trk, const rb::Cluster &slice){
905  //prepare an cluster object to save the shower hits
906  //rb::Cluster shower;
907 
908  //bool isContained = true;
909  //bool isShower = false;
910  bool hasShower = false;
911 
912  //int recoDecayPlane;
913  //recoDecayPlane = decayPlane;
914 
915  //loop over slice cell hits
916  //if(j>=slices.size()){
917  //continue;
918  //}
919 
920  //get the direction of this track
921  TVector3 direction = trk->Dir();
922 
923  //double length = trk->TotalLength();
924 
925  int minPlane = trk->MinPlane();
926  int maxPlane = trk->MaxPlane();
927  int nPlane = geom->NPlanes();
928 
929  if(direction.Z()<0.75 ||
930  fabs(maxPlane-minPlane)<0.3*nPlane ||
931  hasShower = false)//if this is simply not the track we want, remove all the hits in the slice
932 
933  for(unsigned int iCell=0; iCell<slice->NCell(); ++iCell){
934  art::Ptr<rb::CellHit> chit = slices[j]->Cell(iCell);
935  const rb::RecoHit rhit = slices[j]->RecoHit(slices[j]->Cell(iCell));
936 
937  //edepSlice[chit->Plane()][chit->Cell()] += rhit.GeV();
938  //need to know where is this cell hit
939  double xyz[3] = {0.};
940  geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
941 
942  //use truth for now
943  //isShower = (particle->Momentum().Pz()>0)==(xyz[2]>particle->Position().Z());
944  isShower = (trk->Dir().Z()>0)==(chit->Plane()>recoDecayPlane);
945  if(!isShower) continue;
946  hasShower = true;
947 
948  isContained*=geom->isInsideFiducialVolume(xyz[0], xyz[1], xyz[2]);
949 
950  if(!isContained){
951  break;
952  }
953  edepShower[chit->Plane()][chit->Cell()] += rhit.GeV();
954  shower.Add(chit);
955 
956  }//end loop over cell hits
957 
958  if(!hasShower) continue;
959  if(!isContained) continue;
960 
961  //loop over track cell hits
962  for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){
963  art::Ptr<rb::CellHit> chit = trk->Cell(iCell);
964  const rb::RecoHit rhit = trk->RecoHit(trk->Cell(iCell));
965 
966  edepTrack[chit->Plane()][chit->Cell()] += rhit.GeV();
967  }//end loop over cell hits
968 
969  //use true vertex for now
970  TVector3 eStart(particle->Position().X(), particle->Position().Y(), particle->Position().Z());
971 
972  rb::Track etrack(shower, eStart, trk->Dir());
973  etrack.AppendTrajectoryPoint(trk->Stop());
974 
975  showercol->push_back(shower);
976  etrackcol->push_back(etrack);
977  dokeep = true;
978  fTree->Fill();
979  break;
980  }//end loop over tracks
981 }
982 */
983 
984 //find the DiF shower region according to the truth
986 
988  const sim::ParticleNavigator &pnav = bt->ParticleNavigator();
990 
991  if(!(bt->HaveTruthInfo())){
992  mf::LogWarning("BremShowerFilter") <<"No truth info found! ";
993  return;
994  }
995 
996  //int nPlane = geom->NPlanes();
997 
998  if(fabs(trk->Dir().Z())<fAngleCut) return;
999  if(fabs(trk->MaxPlane()-trk->MinPlane())<fPlaneCut) return;
1000 
1001  //loop over true particles to find DiF
1002  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i) {
1003  const sim::Particle *particle = (*i).second;
1004 
1005  //find what is the ture particle of the reco track by timing
1006  //if(particle->Position().T()>trk->MinTNS() &&
1007  // particle->Position().T()<trk->MaxTNS() ){//found it!
1008  //}
1009  //else continue;//don't look at the rest!
1010 
1011  if(particle->Position().T()<trk->MinTNS() ||
1012  particle->Position().T()>trk->MaxTNS() )
1013  continue;
1014 
1015  //----------------------------------------------------------
1016  //is this a DiF
1017  //bool isDiF = false;
1018 
1019  // electrons only
1020  if ( abs(particle->PdgCode()) != 11 )
1021  continue;
1022 
1023  // must from decay
1024  if ( particle->Process() != "Decay" )
1025  continue;
1026 
1027  // make sure electrons start inside the detector
1028  if (!geom->isInsideFiducialVolume(particle->Position().X(),
1029  particle->Position().Y(),
1030  particle->Position().Z()))
1031  continue;
1032 
1033  // must decay from muons
1034  const sim::Particle *mother = pnav[particle->Mother()];
1035  if(abs(mother->PdgCode()) != 13) continue;
1036 
1037  //particleId = particle->Mother();//this is suppose to be the id of muon
1038 
1039 
1040  // not Michel electrons
1041  if ( particle->Momentum().E() < 0.06) continue;
1042 
1043  //get flshits of el and mu
1044  if (bt->ParticleToFLSHit(particle->TrackId()).size()==0) continue;
1045  //if (bt->ParticleToFLSHit(particle->Mother()).size()==0) continue;
1046  //we found DiF!
1047  //isDiF = true;
1048 
1049  if (particle->Momentum().E() < fEnergyCut) continue;
1050 
1051 
1052  /*p4Mu->SetPxPyPzE(mother->Momentum().Px(),
1053  mother->Momentum().Py(),
1054  mother->Momentum().Pz(),
1055  mother->Momentum().E());
1056 
1057  v4->SetXYZT(particle->Position().X(),
1058  particle->Position().Y(),
1059  particle->Position().Z(),
1060  particle->Position().T());
1061 
1062  p4El->SetPxPyPzE(particle->Momentum().Px(),
1063  particle->Momentum().Py(),
1064  particle->Momentum().Pz(),
1065  particle->Momentum().E());
1066  */
1067 
1068  const std::vector<sim::FLSHit>& flshitsEl
1069  = bt->ParticleToFLSHit(particle->TrackId());
1070 
1071  int minPlane = 1000;
1072  int maxPlane = 0;
1073 
1074  for(unsigned int h = 1; h < flshitsEl.size(); ++h){
1075  int plane=flshitsEl[h].GetPlaneID();
1076  if(plane<minPlane) minPlane = plane;
1077  if(plane>maxPlane) maxPlane = plane;
1078  }
1079 
1080  startPlane = minPlane;
1081  endPlane = maxPlane;
1082 
1083  }//end the loop over true particles
1084 
1085  return;
1086 }//end findShowerByTruth
1087 
1088 namespace bsf{
1090 }
size_t NTrajectoryPoints() const
Definition: Track.h:83
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
double HalfD() const
Definition: CellGeo.cxx:205
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
set< int >::iterator it
int32_t TDC() const
The time of the last baseline sample.
Definition: RawDigit.h:94
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
const Var weight
unsigned short Plane() const
Definition: CellHit.h:39
int Mother() const
Definition: MCParticle.h:212
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
list_type::const_iterator const_iterator
Vertical planes which measure X.
Definition: PlaneGeo.h:28
A collection of associated CellHits.
Definition: Cluster.h:47
void abs(TH1 *hist)
bool trackEndContained(art::Ptr< rb::Track > &trk)
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
std::string Process() const
Definition: MCParticle.h:214
Particle class.
rb::CellHit MakeCellHit(const rawdata::RawDigit *rawdigit)
int TrackId() const
Definition: MCParticle.h:209
Float_t Y
Definition: plot.C:38
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
Float_t Z
Definition: plot.C:38
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
double MinTNS() const
Definition: Cluster.cxx:482
bool trackStartContained(art::Ptr< rb::Track > &trk)
uint32_t getNPretriggeredSamples(const version_t ver) const
Get number of pretriggered samples.
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
int evt
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
void findShowerByReco(art::Ptr< rb::Track > &trk, int &startPlane, int &endPlane)
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
std::string fTrackLabel
label of module making fls hits
double DetHalfHeight() const
Perform a "2 point" Hough transform on a collection of hits.
Definition: View.py:1
z
Definition: test.py:28
double MaxTNS() const
Definition: Cluster.cxx:528
Definition: run.py:1
const art::ProductToken< std::vector< rb::Track > > fTrackToken
void reconfigure(const fhicl::ParameterSet &pset)
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
EventNumber_t event() const
Definition: EventID.h:116
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
::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
float GeV() const
Definition: RecoHit.cxx:69
double DetHalfWidth() const
BremShowerFilter(fhicl::ParameterSet const &pset)
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
T const * get() const
Definition: Ptr.h:321
art::Ptr< rb::CellHit > hit
Definition: WeightedHit.h:23
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
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
Simple hit+weight pair, returned from rb::Cluster::WeightedHits.
Definition: WeightedHit.h:18
Double_t ymin
Definition: plot.C:24
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
double showerStart[3]
const art::ProductToken< std::vector< rb::Cluster > > fClusterToken
float PECorr() const
Definition: RecoHit.cxx:47
void findShowerByTruth(art::Ptr< rb::Track > &trk, int &startPlane, int &endPlane)
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
Float_t X
Definition: plot.C:38
virtual void Clear()
Forget about all owned cell hits.
Definition: Cluster.cxx:279
ProductToken< T > consumes(InputTag const &)
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)
bool inFiducial(double x, double y, double z)
bool isInsideFiducialVolume(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Fiducial Volume?
virtual bool filter(art::Event &evt)
const art::ProductToken< std::vector< rb::CellHit > > fCellHitToken