simfilter::FilterGenInTime Class Reference

Inheritance diagram for simfilter::FilterGenInTime:

Public Member Functions
	FilterGenInTime (fhicl::ParameterSet const &pset)
bool	filter (art::Event &)
virtual void	beginJob ()

Private Member Functions
bool	KeepParticle (simb::MCParticle const &part) const

Private Attributes
std::vector< std::array < double, 6 > >	fCryostatBoundaries
	boundaries of each cryostat More...
double	fMinKE
bool	fKeepOnlyMuons
double	fMinT
double	fMaxT
bool	fSortParticles
bool	fAlwaysPass

Detailed Description

Definition at line 38 of file FilterGenInTime_module.cc.

Constructor & Destructor Documentation

simfilter::FilterGenInTime::FilterGenInTime ( fhicl::ParameterSet const &  pset )

explicit

Definition at line 64 of file FilterGenInTime_module.cc.

                                                                :
    EDFilter{pset},
    fMinKE    (pset.get< double > ("MinEnergy"  , 0.0)       )
    , fKeepOnlyMuons    (pset.get< bool > ("KeepOnlyMuons", false)      )
    , fMinT    (pset.get< double > ("MinT",0.0)      )
    , fMaxT    (pset.get< double > ("MaxT")      )
    , fSortParticles ( pset.get< bool > ("SortParticles",false) )
    , fAlwaysPass (pset.get<bool>("AlwaysPass",false))
  {
    if(fSortParticles) {
      produces< std::vector<simb::MCTruth> >("intime");
      produces< std::vector<simb::MCTruth> >("outtime"); 
    }
    std::cout << "New Filter!\n";

  }

Member Function Documentation

void simfilter::FilterGenInTime::beginJob ( )

virtual

Definition at line 81 of file FilterGenInTime_module.cc.

                                {
    auto const& geom = *art::ServiceHandle<geo::Geometry const>();
    for (auto const &cryo: geom.IterateCryostats()) {
      std::array<double, 6> this_cryo_boundaries {};
      cryo.Boundaries(&this_cryo_boundaries[0]);
      fCryostatBoundaries.push_back(this_cryo_boundaries);
    }
  }

bool simfilter::FilterGenInTime::filter

(

art::Event &

evt

)

Definition at line 222 of file FilterGenInTime_module.cc.

                                           {

    std::unique_ptr< std::vector<simb::MCTruth> > truthInTimePtr(new std::vector<simb::MCTruth>(1));
    std::unique_ptr< std::vector<simb::MCTruth> > truthOutOfTimePtr(new std::vector<simb::MCTruth>(1));

    simb::MCTruth & truthInTime = truthInTimePtr->at(0);
    simb::MCTruth & truthOutOfTime = truthOutOfTimePtr->at(0);


    //get the list of particles from this event
    art::ServiceHandle<geo::Geometry const> geom;

    //std::vector< art::Handle< std::vector<simb::MCTruth> > > allmclists;
    //evt.getManyByType(allmclists);
    auto allmclists = evt.getMany< std::vector<simb::MCTruth> >();
    bool keepEvent=false;
    for(size_t mcl = 0; mcl < allmclists.size(); ++mcl){
      art::Handle< std::vector<simb::MCTruth> > mclistHandle = allmclists[mcl];
      for(size_t m = 0; m < mclistHandle->size(); ++m){
        art::Ptr<simb::MCTruth> mct(mclistHandle, m);

        for(int ipart=0;ipart<mct->NParticles();ipart++){
          bool kp=KeepParticle(mct->GetParticle(ipart));

          if(kp
             && (!fKeepOnlyMuons || abs(mct->GetParticle(ipart).PdgCode())==13 )
             && mct->GetParticle(ipart).E()-mct->GetParticle(ipart).Mass()>fMinKE){
            keepEvent = true;
            if(!fSortParticles) break;
          }

          if(fSortParticles){
            simb::MCParticle particle = mct->GetParticle(ipart);
            if(kp) truthInTime.Add(particle);
            if(!kp) truthOutOfTime.Add(particle);
          }

        }//end loop over particles

        if(!fSortParticles && keepEvent) break;

      }//end loop over mctruth col

      if(!fSortParticles && keepEvent) break;

    }//end loop over all mctruth lists

    if(fSortParticles){
      evt.put(std::move(truthInTimePtr),"intime");
      evt.put(std::move(truthOutOfTimePtr),"outtime");
    }

    return (keepEvent || fAlwaysPass);
  }

bool simfilter::FilterGenInTime::KeepParticle ( simb::MCParticle const &  part ) const

private

Definition at line 91 of file FilterGenInTime_module.cc.

                                                                     {
    const TLorentzVector& v4 = part.Position();
    const TLorentzVector& p4 = part.Momentum();
    // origin of particle
    double x0[3] = {v4.X(),  v4.Y(),  v4.Z() };
    // normalized direction of particle
    double dx[3] = {p4.Px() / p4.Vect().Mag(), p4.Py() / p4.Vect().Mag(), p4.Pz() / p4.Vect().Mag()};

    // tolernace for treating number as "zero"
    double eps = 1e-5;

    //check to see if particle crosses boundary of any cryostat within appropriate time window
    std::vector<bool> intersects_cryo(fCryostatBoundaries.size(), false);
    std::vector<bool> inside_cryo(fCryostatBoundaries.size(), false);
    std::vector<double> distance_to_cryo(fCryostatBoundaries.size(), 0.);

    // Check to see if particle intersects any cryostat
    //
    // Algorithmically, this is looking for ray-box intersection. This is a common problem in 
    // computer graphics. The algorithm below is taken from "Graphics Gems", Academic Press, 1990
    for (size_t i_cryo = 0; i_cryo < fCryostatBoundaries.size(); i_cryo++) {
      auto const &bound = fCryostatBoundaries[i_cryo];
      std::array<int, 3> quadrant {}; // 0 == RIGHT, 1 == LEFT, 2 == MIDDLE
      std::array<double, 3> candidatePlane {};
      std::array<double, 3> coord {};

      std::array<double, 3> bound_lo = {{bound[0], bound[2], bound[4]}};
      std::array<double, 3> bound_hi = {{bound[1], bound[3], bound[5]}};

      // First check if origin is inside box
      // Also check which of the two planes in each dimmension is the 
      // "candidate" for the ray to hit
      bool inside = true;
      for (int i = 0; i < 3; i++) {
        if (x0[i] < bound_lo[i]) {
          quadrant[i] = 1; // LEFT
          candidatePlane[i] = bound_lo[i];
          inside = false;
        }
        else if (x0[i] > bound_hi[i]) {
          quadrant[i] = 0; // RIGHT
          candidatePlane[i] = bound_hi[i];
          inside = false;
        }
        else {
          quadrant[i] = 2; // MIDDLE
        }
      }

      if (inside) {
        inside_cryo[i_cryo] = true;
        // if we're inside the cryostat, then we do intersect it
        intersects_cryo[i_cryo] = true;
        continue;
      }

      // ray origin is outside the box -- calculate the distance to the cryostat and see if it intersects
      
      // calculate distances to candidate planes
      std::array<double, 3> maxT {};
      for (int i = 0; i < 3; i++) {
        if (quadrant[i] != 2 /* MIDDLE */ && abs(dx[i]) > eps) {
          maxT[i] = (candidatePlane[i] - x0[i]) / dx[i];
        }
        // if a ray origin is between two the two planes in a dimmension, it would never hit that plane first
        else {
          maxT[i] = -1;
        }
      }

      // The plane on the box that the ray hits is the one with the largest distance
      int whichPlane = 0;
      for (int i = 1; i < 3; i++) {
        if (maxT[whichPlane] < maxT[i]) whichPlane = i;
      }

      // check if the candidate intersection point is inside the box

      // no intersection
      if (maxT[whichPlane] < 0.) {
        intersects_cryo[i_cryo] = false;
        continue;
      }

      for (int i = 0; i < 3; i++) {
        if (whichPlane != i) {
          coord[i] = x0[i] + maxT[whichPlane] * dx[i];
        }
        else {
          coord[i] = candidatePlane[i];
        }
      }


      // check if intersection is in box
      intersects_cryo[i_cryo] = true;
      for (int i = 0; i < 3; i++) {
        if (coord[i] < bound_lo[i] || coord[i] > bound_hi[i]) {
          intersects_cryo[i_cryo] = false;
        }
      }

      if (intersects_cryo[i_cryo]) {
        distance_to_cryo[i_cryo] = maxT[whichPlane];
      }
    }

    // check if any cryostats are intersected in-time
    for (size_t i_cryo = 0; i_cryo < fCryostatBoundaries.size(); i_cryo++) {

      // If the particle originates inside the cryostat, then
      // we can't really say when it will leave. Thus, accept
      // the particle
      if (inside_cryo[i_cryo]) {
        return true;
      }
      // otherwise check arrival time at boundary of cryostat
      if (intersects_cryo[i_cryo]){
        double ptime = (distance_to_cryo[i_cryo] * 1e-2 /* cm -> m */) / (TMath::C()*sqrt(1-pow(part.Mass()/part.E(),2))) /* velocity */;
        double totT=part.T()+ptime*1e9 /* s -> ns */;
        if(totT>fMinT && totT<fMaxT){
          return true;
        }
      }
    }



    return false;
  }

Member Data Documentation

bool simfilter::FilterGenInTime::fAlwaysPass

private

Definition at line 57 of file FilterGenInTime_module.cc.

std::vector<std::array<double, 6> > simfilter::FilterGenInTime::fCryostatBoundaries

private

boundaries of each cryostat

Definition at line 52 of file FilterGenInTime_module.cc.

bool simfilter::FilterGenInTime::fKeepOnlyMuons

private

Definition at line 54 of file FilterGenInTime_module.cc.

double simfilter::FilterGenInTime::fMaxT

private

Definition at line 55 of file FilterGenInTime_module.cc.

double simfilter::FilterGenInTime::fMinKE

private

Definition at line 53 of file FilterGenInTime_module.cc.

double simfilter::FilterGenInTime::fMinT

private

Definition at line 55 of file FilterGenInTime_module.cc.

bool simfilter::FilterGenInTime::fSortParticles

private

Definition at line 56 of file FilterGenInTime_module.cc.

---

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

• FilterGenInTime_module.cc