264 art::Handle<std::vector<simb::MCParticle>> mcParticleHandle;
268 art::Handle<std::vector<simb::MCTruth>> mctruthListHandle;
272 art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
274 art::ServiceHandle<geo::Geometry const> geom;
275 auto const clock_data =
276 art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(
evt);
277 auto const det_prop =
278 art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(
evt, clock_data);
281 art::Handle<std::vector<recob::EndPoint2D>> vertex2dHandle;
285 art::Handle<std::vector<recob::Vertex>> vertex3dListHandle;
289 std::vector<double> WirePitch_CurrentPlane(geom->Views().size(),
294 for (
auto v : geom->Views()) {
295 WirePitch_CurrentPlane[vn] = geom->WirePitch(v);
302 float DriftVelocity = det_prop.DriftVelocity(det_prop.Efield(), det_prop.Temperature());
303 float TimetoCm = TimeTick * DriftVelocity;
308 art::PtrVector<simb::MCTruth> mclist;
309 for (
unsigned int ii = 0; ii < mctruthListHandle->size(); ++ii) {
310 art::Ptr<simb::MCTruth> mctparticle(mctruthListHandle, ii);
311 mclist.push_back(mctparticle);
315 float truth_vertex[5] = {0.};
317 uint32_t VtxWireNum[3] = {0};
318 double VtxTimeTick[3] = {0.};
320 double VtxWireNum_InCM[3] = {0.};
321 double VtxTimeTick_InCM[3] = {0.};
324 for (
unsigned int i = 0; i < mclist.size(); ++i) {
325 art::Ptr<simb::MCTruth> mc(mclist[i]);
326 simb::MCParticle neut(mc->GetParticle(i));
329 truth_vertex[0] = neut.Vx();
330 truth_vertex[1] = neut.Vy();
331 truth_vertex[2] = neut.Vz();
341 for (
auto const& pid : geom->IteratePlaneIDs()) {
344 VtxWireNum[pid.Plane] = geom->NearestWire(truth_vertex, pid.Plane, pid.TPC, pid.Cryostat);
347 mf::LogWarning(
"FeatureVertexFinderAna") <<
"Can't find nearest wire";
350 VtxTimeTick[pid.Plane] =
351 det_prop.ConvertXToTicks(truth_vertex[0], pid.Plane, pid.TPC, pid.Cryostat) +
352 det_prop.GetXTicksOffset(pid.Plane, pid.TPC, pid.Cryostat);
355 VtxWireNum_InCM[pid.Plane] = VtxWireNum[pid.Plane] * WirePitch_CurrentPlane[pid.Plane];
356 VtxTimeTick_InCM[pid.Plane] = VtxTimeTick[pid.Plane] * TimetoCm;
376 art::PtrVector<recob::EndPoint2D> vert2d;
379 double Vertex2d_TimeTick[10000] = {
381 double Vertex2d_Wire[10000] = {0.};
383 double Vertex2d_TimeTick_InCM[10000] = {0.};
384 double Vertex2d_Wire_InCM[10000] = {0.};
387 int n2dVtxPlane0 = 0, n2dVtxPlane1 = 0, n2dVtxPlane2 = 0;
389 bool vertexWstrengthplane0 =
false,
390 vertexWstrengthplane1 =
394 for (
size_t ii = 0; ii < vertex2dHandle->size(); ++ii) {
395 art::Ptr<recob::EndPoint2D>
vertex(vertex2dHandle, ii);
400 if (vert2d.size() > 0) {
403 for (
auto const& pid : geom->IteratePlaneIDs()) {
404 for (
size_t ww = 0; ww < vert2d.size(); ++ww) {
406 if (vert2d[ww]->
WireID().planeID() != pid) {
continue; }
408 Vertex2d_TimeTick[n2dVtx] = vert2d[ww]->DriftTime();
409 Vertex2d_Wire[n2dVtx] = vert2d[ww]->WireID().Wire;
412 Vertex2d_Wire_InCM[n2dVtx] = Vertex2d_Wire[n2dVtx] * WirePitch_CurrentPlane[pid.Plane];
413 Vertex2d_TimeTick_InCM[n2dVtx] = Vertex2d_TimeTick[n2dVtx] * TimetoCm;
417 float RecoCheck_TimeTick = Vertex2d_TimeTick[n2dVtx] - VtxTimeTick[pid.Plane];
418 float RecoCheck_WireNum = Vertex2d_Wire[n2dVtx] - VtxWireNum[pid.Plane];
420 float RecoCheck_TimeInCm = Vertex2d_TimeTick_InCM[n2dVtx] - VtxTimeTick_InCM[pid.Plane];
421 float RecoCheck_WireInCm = Vertex2d_Wire_InCM[n2dVtx] - VtxWireNum_InCM[pid.Plane];
423 if (vert2d[ww]->Strength() > -1) {
424 if (pid.Plane == 0) {
425 vertexWstrengthplane0 =
true;
443 if (pid.Plane == 1) {
444 vertexWstrengthplane1 =
true;
461 if (pid.Plane == 2) {
491 art::PtrVector<recob::Vertex> Vertexlist;
493 double xyz[3] = {0.};
495 for (
unsigned int ii = 0; ii < vertex3dListHandle->size(); ++ii) {
496 art::Ptr<recob::Vertex> vertex3d(vertex3dListHandle, ii);
497 Vertexlist.push_back(vertex3d);
501 if (Vertexlist.size() > 0 && vertexWstrengthplane0 && vertexWstrengthplane1) {
504 for (
unsigned int ww = 0; ww < Vertexlist.size(); ww++) {
505 Vertexlist[ww]->XYZ(xyz);
510 double DeltaX = xyz[0] - truth_vertex[0];
511 double DeltaY = xyz[1] - truth_vertex[1];
512 double DeltaZ = xyz[2] - truth_vertex[2];
514 double DeltaXoverTrueX = DeltaX / truth_vertex[0];
515 double DeltaYoverTrueY = DeltaY / truth_vertex[0];
516 double DeltaZoverTrueZ = DeltaZ / truth_vertex[0];
TH1F * fTwoDTimeInCmPlane2
TH2D * fRecoCheck3dVtxXvsX
TH1F * fTwoDTimeTickPlane2
std::string fEndPoint2dModuleLabel
TH1F * fTwoDTimeInCmPlane0
TH1F * fTwoDTimeInCmPlane1
TH1F * fTruthTimeInCmPlane0
TH1F * fTwoDWireNumberPlane1
TH1F * fTruthWireNumberPlane0
TH1F * fTruthTimeTickPlane2
std::string fVertexModuleLabel
TH1F * fTwoDWireNumberPlane2
TH1F * fRecoCheck2dWireNumPlane1
TH1F * fTruthWireInCmPlane2
TH1F * fTruthWireInCmPlane1
TH1F * fRecoCheck2dTimeInCmPlane2
TH1F * fTwoDStrengthPlane2
TH1F * fTruthTimeInCmPlane1
TH1F * fRecoCheck2dTimeInCmPlane1
TH1F * fRecoCheck2dTimeTickPlane2
std::string fGenieModuleLabel
TH1F * fRecoCheck2dTimeTickPlane0
std::string fLArG4ModuleLabel
TH2D * fRecoCheck3dVtxZvsZ
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
TH1F * fRecoCheck2dWireInCmPlane2
TH1F * fTwoDStrengthPlane0
TH1F * fTwoDStrengthPlane1
TH1F * fTwoDTimeTickPlane0
TH1F * fTwoDTimeTickPlane1
TH1F * fTruthTimeTickPlane0
TH1F * fRecoCheck2dWireNumPlane2
TH1F * fTruthWireNumberPlane2
TH1F * fRecoCheck2dTimeInCmPlane0
TH1F * fRecoCheck2dWireNumPlane0
TH1F * fTruthTimeTickPlane1
TH1F * fRecoCheck2dWireInCmPlane1
TH1F * fRecoCheck2dTimeTickPlane1
TH1F * fTwoDWireInCmPlane1
TH2D * fRecoCheck3dVtxYvsY
TH1F * fRecoCheck2dWireInCmPlane0
TH1F * fTwoDWireInCmPlane0
TH1F * fTruthWireInCmPlane0
TH1F * fTruthWireNumberPlane1
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
TH1F * fTwoDWireInCmPlane2
TH1F * fTruthTimeInCmPlane2
TH1F * fTwoDWireNumberPlane0