|
| 1 | +//////////////////////////////////////////////////////////////////////// |
| 2 | +// Class: CRTVetoProducer |
| 3 | +// Plugin Type: producer |
| 4 | +// File: CRTVetoProducer_module.cc |
| 5 | +// |
| 6 | +// Author: Alex Antonakis (aantonak@ucsb.edu) |
| 7 | +//////////////////////////////////////////////////////////////////////// |
| 8 | + |
| 9 | +#include "art/Framework/Core/EDProducer.h" |
| 10 | +#include "art/Framework/Core/ModuleMacros.h" |
| 11 | +#include "art/Framework/Principal/Event.h" |
| 12 | +#include "art/Framework/Principal/Handle.h" |
| 13 | +#include "art/Framework/Principal/Run.h" |
| 14 | +#include "art/Framework/Principal/SubRun.h" |
| 15 | +#include "canvas/Utilities/InputTag.h" |
| 16 | +#include "canvas/Persistency/Common/FindOneP.h" |
| 17 | +#include "fhiclcpp/ParameterSet.h" |
| 18 | +#include "messagefacility/MessageLogger/MessageLogger.h" |
| 19 | + |
| 20 | +#include "Math/GenVector/PositionVector2D.h" |
| 21 | +#include "Math/GenVector/DisplacementVector2D.h" |
| 22 | + |
| 23 | +#include "lardata/Utilities/AssociationUtil.h" |
| 24 | + |
| 25 | +#include "sbnobj/SBND/CRT/CRTCluster.hh" |
| 26 | +#include "sbnobj/SBND/CRT/CRTSpacePoint.hh" |
| 27 | +#include "sbnobj/SBND/CRT/CRTTrack.hh" |
| 28 | +#include "sbnobj/SBND/CRT/CRTVeto.hh" |
| 29 | + |
| 30 | +#include "sbndcode/Geometry/GeometryWrappers/CRTGeoAlg.h" |
| 31 | +#include "sbndcode/CRT/CRTUtils/CRTCommonUtils.h" |
| 32 | + |
| 33 | + |
| 34 | +#include "TMath.h" |
| 35 | + |
| 36 | +#include <memory> |
| 37 | + |
| 38 | + |
| 39 | +namespace sbnd::crt { |
| 40 | + class CRTVetoProducer; |
| 41 | +} |
| 42 | + |
| 43 | + |
| 44 | +class sbnd::crt::CRTVetoProducer : public art::EDProducer { |
| 45 | +public: |
| 46 | + explicit CRTVetoProducer(fhicl::ParameterSet const& p); |
| 47 | + |
| 48 | + CRTVetoProducer(CRTVetoProducer const&) = delete; |
| 49 | + CRTVetoProducer(CRTVetoProducer&&) = delete; |
| 50 | + CRTVetoProducer& operator=(CRTVetoProducer const&) = delete; |
| 51 | + CRTVetoProducer& operator=(CRTVetoProducer&&) = delete; |
| 52 | + |
| 53 | + void produce(art::Event& e) override; |
| 54 | + |
| 55 | + // Functions to check if a point is in a specific Tagger |
| 56 | + bool inSouth(const double &x, const double &y, const double &z); |
| 57 | + bool inNorth(const double &x, const double &y, const double &z); |
| 58 | + bool inEast(const double &x, const double &y, const double &z); |
| 59 | + bool inWest(const double &x, const double &y, const double &z); |
| 60 | + bool inTopLow(const double &x, const double &y, const double &z); |
| 61 | + bool inTopHigh(const double &x, const double &y, const double &z); |
| 62 | + |
| 63 | +private: |
| 64 | + |
| 65 | + CRTGeoAlg fCRTGeoAlg; |
| 66 | + double fWindowStart; |
| 67 | + double fWindowEnd; |
| 68 | + std::string fCRTSpacePointModuleLabel; |
| 69 | +}; |
| 70 | + |
| 71 | + |
| 72 | +sbnd::crt::CRTVetoProducer::CRTVetoProducer(fhicl::ParameterSet const& p) |
| 73 | + : EDProducer{p} |
| 74 | + , fCRTGeoAlg(p.get<fhicl::ParameterSet>("CRTGeoAlg")) |
| 75 | + , fWindowStart(p.get<double>("WindowStart")) |
| 76 | + , fWindowEnd(p.get<double>("WindowEnd")) |
| 77 | + , fCRTSpacePointModuleLabel(p.get<std::string>("CRTSpacePointModuleLabel")) |
| 78 | +{ |
| 79 | + produces<CRTVeto>(); |
| 80 | + // produces<art::Assns<CRTSpacePoint, CRTVeto>>(); |
| 81 | +} |
| 82 | + |
| 83 | +void sbnd::crt::CRTVetoProducer::produce(art::Event& e) |
| 84 | +{ |
| 85 | + //auto crtveto = std::make_unique<CRTVeto>(); |
| 86 | + //auto vetoSpacePointAssn = std::make_unique<art::Assns<CRTSpacePoint, CRTVeto>>(); |
| 87 | + |
| 88 | + // Get the Space Points for this event |
| 89 | + art::Handle<std::vector<CRTSpacePoint>> CRTSpacePointHandle; |
| 90 | + e.getByLabel(fCRTSpacePointModuleLabel, CRTSpacePointHandle); |
| 91 | + if(!CRTSpacePointHandle.isValid()){ |
| 92 | + std::cout << "CRTSpacePoint product " << fCRTSpacePointModuleLabel << " not found..." << std::endl; |
| 93 | + throw std::exception(); |
| 94 | + } |
| 95 | + std::vector<art::Ptr<CRTSpacePoint>> CRTSpacePointVec; |
| 96 | + art::fill_ptr_vector(CRTSpacePointVec, CRTSpacePointHandle); |
| 97 | + |
| 98 | + int S = 0; |
| 99 | + int N = 0; |
| 100 | + int E = 0; |
| 101 | + int W = 0; |
| 102 | + int TL = 0; |
| 103 | + int TH = 0; |
| 104 | + |
| 105 | + // loop over the space points |
| 106 | + for (unsigned int i = 0; i < CRTSpacePointVec.size(); i++){ |
| 107 | + |
| 108 | + double x = CRTSpacePointVec[i]->X(); |
| 109 | + double y = CRTSpacePointVec[i]->Y(); |
| 110 | + double z = CRTSpacePointVec[i]->Z(); |
| 111 | + double t = CRTSpacePointVec[i]->Ts0(); |
| 112 | + // check if in the beam window --> Make configureable |
| 113 | + if ((t > fWindowEnd) || (t < fWindowStart)) { |
| 114 | + continue; |
| 115 | + } |
| 116 | + // We have a valid space point --> make assn |
| 117 | + //util::CreateAssn(*this, e, *crtveto, CRTSpacePointVec[i], *vetoSpacePointAssn); |
| 118 | + |
| 119 | + if (inSouth(x, y, z)) { |
| 120 | + S += 1; |
| 121 | + } |
| 122 | + else if (inNorth(x, y, z)) { |
| 123 | + N += 1; |
| 124 | + } |
| 125 | + else if (inEast(x, y, z)) { |
| 126 | + E += 1; |
| 127 | + } |
| 128 | + else if (inWest(x, y, z)) { |
| 129 | + W += 1; |
| 130 | + } |
| 131 | + else if (inTopLow(x, y, z)) { |
| 132 | + TL += 1; |
| 133 | + } |
| 134 | + else if (inTopHigh(x, y, z)) { |
| 135 | + TH += 1; |
| 136 | + } |
| 137 | + else { |
| 138 | + std::cout << "Not a valid SpacePoint --> Outside Geometry !!!" << std::endl; |
| 139 | + } |
| 140 | + |
| 141 | + } // end loop over space points |
| 142 | + |
| 143 | + int Nwalls_all = S + N + E + W; |
| 144 | + int Nwalls_noNorth = S + E + W; |
| 145 | + int Ntop = TL + TH; |
| 146 | + |
| 147 | + bool v0 = false; |
| 148 | + bool v1 = false; |
| 149 | + bool v2 = false; |
| 150 | + bool v3 = false; |
| 151 | + |
| 152 | + // V0 |
| 153 | + if ((Nwalls_all + Ntop) > 0) { |
| 154 | + v0 = true; |
| 155 | + } |
| 156 | + // V1 |
| 157 | + if ( (Nwalls_all > 0) || ((TL > 0) && (TH > 0)) ) { |
| 158 | + v1 = true; |
| 159 | + } |
| 160 | + // V2 |
| 161 | + if ((Nwalls_noNorth + Ntop) > 0) { |
| 162 | + v2 = true; |
| 163 | + } |
| 164 | + // V3 |
| 165 | + if ( (Nwalls_noNorth > 0) || ((TL > 0) && (TH > 0)) ) { |
| 166 | + v3 = true; |
| 167 | + } |
| 168 | + |
| 169 | + //crtveto = CRTVeto(v0, v1, v2, v3); |
| 170 | + auto crtveto = std::make_unique<CRTVeto>(v0, v1, v2, v3); |
| 171 | + |
| 172 | + e.put(std::move(crtveto)); |
| 173 | + //e.put(std::move(vetoSpacePointAssn)); |
| 174 | + |
| 175 | +} // end of produce |
| 176 | + |
| 177 | +bool sbnd::crt::CRTVetoProducer::inSouth(const double &x, const double &y, const double &z) |
| 178 | +{ |
| 179 | + if ((z > -1000) && (z < -179) && (y > -370) && (y < 400)) { |
| 180 | + return true; |
| 181 | + } |
| 182 | + return false; |
| 183 | +} |
| 184 | +bool sbnd::crt::CRTVetoProducer::inNorth(const double &x, const double &y, const double &z) |
| 185 | +{ |
| 186 | + // check that this isn't affected by x from east and west --> TODO |
| 187 | + if ((z > 746) && (y > -370) && (y < 400)) { |
| 188 | + return true; |
| 189 | + } |
| 190 | + return false; |
| 191 | +} |
| 192 | + |
| 193 | +bool sbnd::crt::CRTVetoProducer::inEast(const double &x, const double &y, const double &z) |
| 194 | +{ |
| 195 | + //east_sel = "x < -380 and z < 770 and -370 < y < 400" |
| 196 | + if ((x < -380) && (z < 770) && (y > -370) && (y < 400)) { |
| 197 | + return true; |
| 198 | + } |
| 199 | + return false; |
| 200 | +} |
| 201 | + |
| 202 | +bool sbnd::crt::CRTVetoProducer::inWest(const double &x, const double &y, const double &z) |
| 203 | +{ |
| 204 | + //west_sel = "x < -380 and z < 770 and -370 < y < 400" |
| 205 | + if ((x > 380) && (z < 770) && (y > -370) && (y < 400)) { |
| 206 | + return true; |
| 207 | + } |
| 208 | + return false; |
| 209 | +} |
| 210 | + |
| 211 | +bool sbnd::crt::CRTVetoProducer::inTopLow(const double &x, const double &y, const double &z) |
| 212 | +{ |
| 213 | + if ((y > 400) && (y < 700)) { |
| 214 | + return true; |
| 215 | + } |
| 216 | + return false; |
| 217 | +} |
| 218 | + |
| 219 | +bool sbnd::crt::CRTVetoProducer::inTopHigh(const double &x, const double &y, const double &z) |
| 220 | +{ |
| 221 | + if (y > 700) { |
| 222 | + return true; |
| 223 | + } |
| 224 | + return false; |
| 225 | +} |
| 226 | + |
| 227 | + |
| 228 | + |
| 229 | +DEFINE_ART_MODULE(sbnd::crt::CRTVetoProducer) |
0 commit comments