Public Types | Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | Private Attributes | List of all members
bsf::BremShowerFilter Class Reference
Inheritance diagram for bsf::BremShowerFilter:
art::EDFilter art::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Public Types

using ModuleType = EDFilter
 
using WorkerType = WorkerT< EDFilter >
 
template<typename UserConfig >
using Table = ProducerBase::Table< UserConfig >
 

Public Member Functions

 BremShowerFilter (fhicl::ParameterSet const &pset)
 
virtual ~BremShowerFilter ()
 
virtual bool filter (art::Event &evt)
 
void reconfigure (const fhicl::ParameterSet &pset)
 
void beginJob ()
 
void findShowerByReco (art::Ptr< rb::Track > &trk, int &startPlane, int &endPlane)
 
void findShowerByTruth (art::Ptr< rb::Track > &trk, int &startPlane, int &endPlane)
 
bool inFiducial (double x, double y, double z)
 
bool trackStartContained (art::Ptr< rb::Track > &trk)
 
bool trackEndContained (art::Ptr< rb::Track > &trk)
 
template<typename PROD , BranchType B = InEvent>
ProductID getProductID (std::string const &instanceName={}) const
 
template<typename PROD , BranchType B>
ProductID getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
 
bool modifiesEvent () const
 
template<typename T , BranchType = InEvent>
ProductToken< Tconsumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< Tconsumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TconsumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< TmayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< TmayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TmayConsumeView (InputTag const &it)
 
base_engine_tcreateEngine (seed_t seed)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumernon_module_context ()
 

Static Public Attributes

static constexpr bool Pass {true}
 
static constexpr bool Fail {false}
 

Protected Member Functions

CurrentProcessingContext const * currentContext () const
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

Private Attributes

bool dokeep
 
const art::ProductToken< std::vector< rb::Track > > fTrackToken
 
const art::ProductToken< std::vector< rb::Cluster > > fClusterToken
 
const art::ProductToken< std::vector< rb::CellHit > > fCellHitToken
 
std::string fTrackLabel
 label of module making fls hits More...
 
bool fWriteParentSliceMode
 
bool fApplyFilter
 
bool fShowerByTruth
 
bool fShowerByReco
 
bool fFindDiF
 
float fAngleCut
 
float fPlaneCut
 
float fEnergyCut
 
float fFiducialX
 
float fFiducialY
 
float fFiducialZ
 
int fPlaneMin
 
int fPlaneMax
 
TTree * fTree
 
int run
 
int subrun
 
int event
 
int particleId
 
int trackId
 
int startPlane
 
int endPlane
 
double cosz
 
double eshower
 
double planeDedx [2]
 
double planeDedxMR [2]
 
double cellDedx [2]
 
double cellDedxMR [2]
 
int fDCSDistance
 
TVector3 * vtx = new TVector3(0,0,0)
 
TH1F * fHNTrk
 
TH1F * fHCosZ
 
TH1F * fHCosZCut [5]
 
TH1F * fHStartContainment
 
TH1F * fHEndContainment
 
TH1F * fHNPlane
 
TH1F * fHEShower
 
TH1F * fHTrackStartZ
 
TH1F * fHEdep
 
TH1F * fHADC
 
TH1F * fHPE
 
TH1F * fHPECorr
 
TH2F * fHHitMap
 
TH2F * fHHitMapRW
 

Detailed Description

Definition at line 47 of file BremShowerFilter_module.cc.

Member Typedef Documentation

using art::EDFilter::ModuleType = EDFilter
inherited

Definition at line 37 of file EDFilter.h.

template<typename UserConfig >
using art::EDFilter::Table = ProducerBase::Table<UserConfig>
inherited

Definition at line 46 of file EDFilter.h.

using art::EDFilter::WorkerType = WorkerT<EDFilter>
inherited

Definition at line 38 of file EDFilter.h.

Constructor & Destructor Documentation

bsf::BremShowerFilter::BremShowerFilter ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 155 of file BremShowerFilter_module.cc.

References reconfigure().

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  }
const art::ProductToken< std::vector< rb::Track > > fTrackToken
void reconfigure(const fhicl::ParameterSet &pset)
const art::ProductToken< std::vector< rb::Cluster > > fClusterToken
ProductToken< T > consumes(InputTag const &)
const art::ProductToken< std::vector< rb::CellHit > > fCellHitToken
enum BeamMode string
bsf::BremShowerFilter::~BremShowerFilter ( )
virtual

Definition at line 171 of file BremShowerFilter_module.cc.

171 { }

Member Function Documentation

void bsf::BremShowerFilter::beginJob ( )
virtual

Reimplemented from art::EDFilter.

Definition at line 176 of file BremShowerFilter_module.cc.

References cellDedx, cellDedxMR, cosz, endPlane, eshower, event, fHADC, fHCosZ, fHCosZCut, fHEdep, fHEndContainment, fHEShower, fHHitMap, fHHitMapRW, fHNPlane, fHNTrk, fHPE, fHPECorr, fHStartContainment, fHTrackStartZ, fTree, MECModelEnuComparisons::i, art::TFileDirectory::make(), particleId, planeDedx, planeDedxMR, startPlane, subrun, trackId, and vtx.

176  {
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  }
Definition: run.py:1
T * make(ARGS...args) const
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::consumes ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::consumes ( InputTag const &  it)
inherited

Definition at line 146 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

147 {
148  if (!moduleContext_)
149  return ProductToken<T>::invalid();
150 
151  consumables_[BT].emplace_back(ConsumableType::Product,
152  TypeID{typeid(T)},
153  it.label(),
154  it.instance(),
155  it.process());
156  return ProductToken<T>{it};
157 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 161 of file Consumer.h.

References T.

162 {
163  if (!moduleContext_)
164  return;
165 
166  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
167 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 171 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

172 {
173  if (!moduleContext_)
174  return ViewToken<T>::invalid();
175 
176  consumables_[BT].emplace_back(ConsumableType::ViewElement,
177  TypeID{typeid(T)},
178  it.label(),
179  it.instance(),
180  it.process());
181  return ViewToken<T>{it};
182 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited
CurrentProcessingContext const* art::EDFilter::currentContext ( ) const
protectedinherited
bool bsf::BremShowerFilter::filter ( art::Event evt)
virtual

Implements art::EDFilter.

Definition at line 255 of file BremShowerFilter_module.cc.

References rawdata::RawDigit::ADC(), rb::Cluster::Add(), rb::CellHit::Cell(), geo::PlaneGeo::Cell(), rb::Cluster::Cell(), rb::Cluster::Clear(), cosz, dokeep, endPlane, art::EventID::event(), stan::math::fabs(), fAngleCut, fApplyFilter, fCellHitToken, fClusterToken, fDCSDistance, fHCosZ, fHCosZCut, fHEndContainment, fHNPlane, fHNTrk, fHStartContainment, fHTrackStartZ, findShowerByReco(), findShowerByTruth(), PandAna.Demos.pi0_spectra::fmt, fPlaneCut, fPlaneMax, fPlaneMin, fShowerByReco, fShowerByTruth, fTrackLabel, fTrackToken, fTree, fWriteParentSliceMode, geom(), art::Ptr< T >::get(), art::DataViewImpl::getByToken(), geo::CellGeo::GetCenter(), daqdataformats::NanoSliceVersionConvention::getNPretriggeredSamples(), rb::RecoHit::GeV(), make_syst_table_plots::h, geo::CellGeo::HalfD(), rb::WeightedHit::hit, art::Event::id(), makeTrainCVSamples::int, rb::RecoHit::IsCalibrated(), rb::Cluster::IsNoise(), it, geo::kX, calib::Calibrator::MakeCellHit(), makeBrightnessMap::maxPlane, rb::Cluster::MeanTNS(), rb::Cluster::NCell(), rb::CellHit::Plane(), geo::GeometryBase::Plane(), NDAPDHVSetting::plane, art::Event::put(), rb::RecoHit::RecoHit(), art::Event::run(), PandAna.Demos.tute_pid_validation::slc, startPlane, art::Event::subRun(), subrun, rawdata::RawDigit::TDC(), trackEndContained(), trackId, trackStartContained(), vtx, ana::weight, X, Y, and Z.

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
SubRunNumber_t subRun() const
Definition: Event.h:72
double HalfD() const
Definition: CellGeo.cxx:205
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
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 Var weight
unsigned short Plane() const
Definition: CellHit.h:39
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
Vertical planes which measure X.
Definition: PlaneGeo.h:28
A collection of associated CellHits.
Definition: Cluster.h:47
bool trackEndContained(art::Ptr< rb::Track > &trk)
const PlaneGeo * Plane(unsigned int i) const
rb::CellHit MakeCellHit(const rawdata::RawDigit *rawdigit)
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
bool trackStartContained(art::Ptr< rb::Track > &trk)
uint32_t getNPretriggeredSamples(const version_t ver) const
Get number of pretriggered samples.
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
Definition: View.py:1
Definition: run.py:1
const art::ProductToken< std::vector< rb::Track > > fTrackToken
EventNumber_t event() const
Definition: EventID.h:116
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
float GeV() const
Definition: RecoHit.cxx:69
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
art::Ptr< rb::CellHit > hit
Definition: WeightedHit.h:23
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
Simple hit+weight pair, returned from rb::Cluster::WeightedHits.
Definition: WeightedHit.h:18
const art::ProductToken< std::vector< rb::Cluster > > fClusterToken
void findShowerByTruth(art::Ptr< rb::Track > &trk, int &startPlane, int &endPlane)
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
EventID id() const
Definition: Event.h:56
const art::ProductToken< std::vector< rb::CellHit > > fCellHitToken
void bsf::BremShowerFilter::findShowerByReco ( art::Ptr< rb::Track > &  trk,
int startPlane,
int endPlane 
)

Definition at line 775 of file BremShowerFilter_module.cc.

References rawdata::RawDigit::ADC(), geo::PlaneGeo::Cell(), rb::Cluster::Cell(), rb::Prong::Dir(), endPlane, eshower, stan::math::fabs(), fEnergyCut, fHADC, fHEdep, fHEShower, fHPE, fHPECorr, fPlaneMax, fPlaneMin, geom(), art::Ptr< T >::get(), rb::RecoHit::GeV(), geo::CellGeo::HalfD(), rb::RecoHit::IsCalibrated(), makeBrightnessMap::maxPlane, rb::Cluster::MaxPlane(), rb::Cluster::MinPlane(), rb::Cluster::NCell(), rb::CellHit::PE(), rb::RecoHit::PECorr(), rb::CellHit::Plane(), geo::GeometryBase::Plane(), rb::Cluster::RecoHit(), and slidt::showerStart.

Referenced by filter().

775  {
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...
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
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
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
unsigned short Plane() const
Definition: CellHit.h:39
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const PlaneGeo * Plane(unsigned int i) const
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
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
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
float GeV() const
Definition: RecoHit.cxx:69
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
T const * get() const
Definition: Ptr.h:321
void geom(int which=0)
Definition: geom.C:163
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
double showerStart[3]
float PECorr() const
Definition: RecoHit.cxx:47
void bsf::BremShowerFilter::findShowerByTruth ( art::Ptr< rb::Track > &  trk,
int startPlane,
int endPlane 
)

Definition at line 985 of file BremShowerFilter_module.cc.

References abs(), sim::ParticleNavigator::begin(), DEFINE_ART_MODULE(), rb::Prong::Dir(), sim::ParticleNavigator::end(), stan::math::fabs(), fAngleCut, fEnergyCut, fPlaneCut, geom(), make_syst_table_plots::h, cheat::BackTracker::HaveTruthInfo(), MECModelEnuComparisons::i, geo::GeometryBase::isInsideFiducialVolume(), makeBrightnessMap::maxPlane, rb::Cluster::MaxPlane(), rb::Cluster::MaxTNS(), rb::Cluster::MinPlane(), rb::Cluster::MinTNS(), simb::MCParticle::Momentum(), simb::MCParticle::Mother(), cheat::BackTracker::ParticleNavigator(), cheat::BackTracker::ParticleToFLSHit(), simb::MCParticle::PdgCode(), NDAPDHVSetting::plane, simb::MCParticle::Position(), simb::MCParticle::Process(), gen_flatrecord::size, and simb::MCParticle::TrackId().

Referenced by filter().

985  {
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
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
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
int Mother() const
Definition: MCParticle.h:212
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
void abs(TH1 *hist)
std::string Process() const
Definition: MCParticle.h:214
int TrackId() const
Definition: MCParticle.h:209
double MinTNS() const
Definition: Cluster.cxx:482
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
double MaxTNS() const
Definition: Cluster.cxx:528
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
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
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
bool isInsideFiducialVolume(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Fiducial Volume?
seed_t art::EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited
template<typename PROD , BranchType B>
ProductID art::EDFilter::getProductID ( std::string const &  instanceName = {}) const
inlineinherited

Definition at line 131 of file EDFilter.h.

References art::EDFilter::moduleDescription_.

Referenced by novaddt::HoughTrackMaker::create_associations().

132  {
133  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
134  instanceName);
135  }
ModuleDescription moduleDescription_
Definition: EDFilter.h:124
template<typename PROD , BranchType B>
ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription,
std::string const &  instanceName 
) const
inherited

Definition at line 56 of file ProducerBase.h.

References art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
bool bsf::BremShowerFilter::inFiducial ( double  x,
double  y,
double  z 
)

Definition at line 700 of file BremShowerFilter_module.cc.

References geo::PlaneGeo::Cell(), geo::GeometryBase::DetHalfHeight(), geo::GeometryBase::DetHalfWidth(), fFiducialX, fFiducialY, fFiducialZ, fPlaneMax, fPlaneMin, geom(), geo::CellGeo::GetCenter(), geo::CellGeo::HalfD(), and geo::GeometryBase::Plane().

Referenced by trackEndContained(), and trackStartContained().

700  {
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 }
double HalfD() const
Definition: CellGeo.cxx:205
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const PlaneGeo * Plane(unsigned int i) const
double DetHalfHeight() const
z
Definition: test.py:28
double DetHalfWidth() const
void geom(int which=0)
Definition: geom.C:163
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::mayConsume ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::mayConsume ( InputTag const &  it)
inherited

Definition at line 189 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

190 {
191  if (!moduleContext_)
192  return ProductToken<T>::invalid();
193 
194  consumables_[BT].emplace_back(ConsumableType::Product,
195  TypeID{typeid(T)},
196  it.label(),
197  it.instance(),
198  it.process());
199  return ProductToken<T>{it};
200 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 204 of file Consumer.h.

References T.

205 {
206  if (!moduleContext_)
207  return;
208 
209  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
210 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it)
inherited

Definition at line 214 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

215 {
216  if (!moduleContext_)
217  return ViewToken<T>::invalid();
218 
219  consumables_[BT].emplace_back(ConsumableType::ViewElement,
220  TypeID{typeid(T)},
221  it.label(),
222  it.instance(),
223  it.process());
224  return ViewToken<T>{it};
225 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
bool art::ProducerBase::modifiesEvent ( ) const
inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID(), and string.

41  {
42  return true;
43  }
static cet::exempt_ptr<Consumer> art::Consumer::non_module_context ( )
staticinherited
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited
void bsf::BremShowerFilter::reconfigure ( const fhicl::ParameterSet pset)

Definition at line 235 of file BremShowerFilter_module.cc.

References fAngleCut, fApplyFilter, fDCSDistance, fEnergyCut, fFiducialX, fFiducialY, fFiducialZ, fFindDiF, fPlaneCut, fPlaneMax, fPlaneMin, fShowerByReco, fShowerByTruth, fTrackLabel, fWriteParentSliceMode, fhicl::ParameterSet::get(), and string.

Referenced by BremShowerFilter().

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  }
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::string fTrackLabel
label of module making fls hits
enum BeamMode string
void art::Consumer::showMissingConsumes ( ) const
protectedinherited

Referenced by art::RootOutput::endJob().

bool bsf::BremShowerFilter::trackEndContained ( art::Ptr< rb::Track > &  trk)

Definition at line 734 of file BremShowerFilter_module.cc.

References inFiducial(), rb::Track::NTrajectoryPoints(), rb::Track::TrajectoryPoint(), and ymin.

Referenced by filter().

734  {
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 }
size_t NTrajectoryPoints() const
Definition: Track.h:83
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
Double_t ymin
Definition: plot.C:24
bool inFiducial(double x, double y, double z)
bool bsf::BremShowerFilter::trackStartContained ( art::Ptr< rb::Track > &  trk)

Definition at line 727 of file BremShowerFilter_module.cc.

References inFiducial(), and rb::Prong::Start().

Referenced by filter().

727  {
728 
729  bool isContained = inFiducial(trk->Start().X(), trk->Start().Y(), trk->Start().Z());
730 
731  return isContained;
732 }
virtual TVector3 Start() const
Definition: Prong.h:73
bool inFiducial(double x, double y, double z)
void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited

Member Data Documentation

double bsf::BremShowerFilter::cellDedx[2]
private

Definition at line 105 of file BremShowerFilter_module.cc.

Referenced by beginJob().

double bsf::BremShowerFilter::cellDedxMR[2]
private

Definition at line 106 of file BremShowerFilter_module.cc.

Referenced by beginJob().

double bsf::BremShowerFilter::cosz
private

Definition at line 101 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and filter().

bool bsf::BremShowerFilter::dokeep
private

Definition at line 63 of file BremShowerFilter_module.cc.

Referenced by filter().

int bsf::BremShowerFilter::endPlane
private

Definition at line 99 of file BremShowerFilter_module.cc.

Referenced by beginJob(), filter(), and findShowerByReco().

double bsf::BremShowerFilter::eshower
private

Definition at line 102 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and findShowerByReco().

int bsf::BremShowerFilter::event
private

Definition at line 93 of file BremShowerFilter_module.cc.

Referenced by beginJob().

constexpr bool art::EDFilter::Fail {false}
staticinherited

Definition at line 33 of file EDFilter.h.

Referenced by evgen::GENIEFilter::filter().

float bsf::BremShowerFilter::fAngleCut
private

Definition at line 79 of file BremShowerFilter_module.cc.

Referenced by filter(), findShowerByTruth(), and reconfigure().

bool bsf::BremShowerFilter::fApplyFilter
private

Definition at line 73 of file BremShowerFilter_module.cc.

Referenced by filter(), and reconfigure().

const art::ProductToken<std::vector<rb::CellHit> > bsf::BremShowerFilter::fCellHitToken
private

Definition at line 68 of file BremShowerFilter_module.cc.

Referenced by filter().

const art::ProductToken<std::vector<rb::Cluster> > bsf::BremShowerFilter::fClusterToken
private

Definition at line 67 of file BremShowerFilter_module.cc.

Referenced by filter().

int bsf::BremShowerFilter::fDCSDistance
private

Definition at line 107 of file BremShowerFilter_module.cc.

Referenced by filter(), and reconfigure().

float bsf::BremShowerFilter::fEnergyCut
private

Definition at line 81 of file BremShowerFilter_module.cc.

Referenced by findShowerByReco(), findShowerByTruth(), and reconfigure().

float bsf::BremShowerFilter::fFiducialX
private

Definition at line 82 of file BremShowerFilter_module.cc.

Referenced by inFiducial(), and reconfigure().

float bsf::BremShowerFilter::fFiducialY
private

Definition at line 83 of file BremShowerFilter_module.cc.

Referenced by inFiducial(), and reconfigure().

float bsf::BremShowerFilter::fFiducialZ
private

Definition at line 84 of file BremShowerFilter_module.cc.

Referenced by inFiducial(), and reconfigure().

bool bsf::BremShowerFilter::fFindDiF
private

Definition at line 78 of file BremShowerFilter_module.cc.

Referenced by reconfigure().

TH1F* bsf::BremShowerFilter::fHADC
private

Definition at line 142 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and findShowerByReco().

TH1F* bsf::BremShowerFilter::fHCosZ
private

Definition at line 131 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and filter().

TH1F* bsf::BremShowerFilter::fHCosZCut[5]
private

Definition at line 132 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and filter().

TH1F* bsf::BremShowerFilter::fHEdep
private

Definition at line 141 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and findShowerByReco().

TH1F* bsf::BremShowerFilter::fHEndContainment
private

Definition at line 134 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and filter().

TH1F* bsf::BremShowerFilter::fHEShower
private

Definition at line 136 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and findShowerByReco().

TH2F* bsf::BremShowerFilter::fHHitMap
private

Definition at line 146 of file BremShowerFilter_module.cc.

Referenced by beginJob().

TH2F* bsf::BremShowerFilter::fHHitMapRW
private

Definition at line 147 of file BremShowerFilter_module.cc.

Referenced by beginJob().

TH1F* bsf::BremShowerFilter::fHNPlane
private

Definition at line 135 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and filter().

TH1F* bsf::BremShowerFilter::fHNTrk
private

Definition at line 130 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and filter().

TH1F* bsf::BremShowerFilter::fHPE
private

Definition at line 143 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and findShowerByReco().

TH1F* bsf::BremShowerFilter::fHPECorr
private

Definition at line 144 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and findShowerByReco().

TH1F* bsf::BremShowerFilter::fHStartContainment
private

Definition at line 133 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and filter().

TH1F* bsf::BremShowerFilter::fHTrackStartZ
private

Definition at line 139 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and filter().

float bsf::BremShowerFilter::fPlaneCut
private

Definition at line 80 of file BremShowerFilter_module.cc.

Referenced by filter(), findShowerByTruth(), and reconfigure().

int bsf::BremShowerFilter::fPlaneMax
private

Definition at line 86 of file BremShowerFilter_module.cc.

Referenced by filter(), findShowerByReco(), inFiducial(), and reconfigure().

int bsf::BremShowerFilter::fPlaneMin
private

Definition at line 85 of file BremShowerFilter_module.cc.

Referenced by filter(), findShowerByReco(), inFiducial(), and reconfigure().

bool bsf::BremShowerFilter::fShowerByReco
private

Definition at line 77 of file BremShowerFilter_module.cc.

Referenced by filter(), and reconfigure().

bool bsf::BremShowerFilter::fShowerByTruth
private

Definition at line 76 of file BremShowerFilter_module.cc.

Referenced by filter(), and reconfigure().

std::string bsf::BremShowerFilter::fTrackLabel
private

label of module making fls hits

Definition at line 70 of file BremShowerFilter_module.cc.

Referenced by filter(), and reconfigure().

const art::ProductToken<std::vector<rb::Track> > bsf::BremShowerFilter::fTrackToken
private

Definition at line 66 of file BremShowerFilter_module.cc.

Referenced by filter().

TTree* bsf::BremShowerFilter::fTree
private

Definition at line 89 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and filter().

bool bsf::BremShowerFilter::fWriteParentSliceMode
private

Definition at line 72 of file BremShowerFilter_module.cc.

Referenced by filter(), and reconfigure().

int bsf::BremShowerFilter::particleId
private

Definition at line 94 of file BremShowerFilter_module.cc.

Referenced by beginJob().

constexpr bool art::EDFilter::Pass {true}
staticinherited

Definition at line 32 of file EDFilter.h.

Referenced by evgen::GENIEFilter::filter().

double bsf::BremShowerFilter::planeDedx[2]
private

Definition at line 103 of file BremShowerFilter_module.cc.

Referenced by beginJob().

double bsf::BremShowerFilter::planeDedxMR[2]
private

Definition at line 104 of file BremShowerFilter_module.cc.

Referenced by beginJob().

int bsf::BremShowerFilter::run
private
int bsf::BremShowerFilter::startPlane
private

Definition at line 98 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and filter().

int bsf::BremShowerFilter::subrun
private

Definition at line 92 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and filter().

int bsf::BremShowerFilter::trackId
private

Definition at line 95 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and filter().

TVector3* bsf::BremShowerFilter::vtx = new TVector3(0,0,0)
private

Definition at line 124 of file BremShowerFilter_module.cc.

Referenced by beginJob(), and filter().


The documentation for this class was generated from the following file: