diff --git a/sbncode/Calibration/TrackCaloSkimmer.h b/sbncode/Calibration/TrackCaloSkimmer.h index f9fb1930e..c69d8589d 100644 --- a/sbncode/Calibration/TrackCaloSkimmer.h +++ b/sbncode/Calibration/TrackCaloSkimmer.h @@ -63,6 +63,8 @@ #include "sbnobj/Common/Calibration/TrackCaloSkimmerObj.h" #include "sbnobj/Common/CRT/CRTHitT0TaggingInfo.hh" #include "sbnobj/Common/CRT/CRTHitT0TaggingTruthInfo.hh" +#include "sbnobj/SBND/CRT/CRTTrack.hh" +#include "sbnobj/SBND/CRT/CRTSpacePoint.hh" #include "ITCSSelectionTool.h" @@ -98,9 +100,12 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer { double t0Pandora; double t0CRTTrack; double t0CRTHit; + double t0CRTSpacePoint; bool hasT0Pandora; bool hasT0CRTTrack; bool hasT0CRTHit; + bool hasT0CRTSpacePoint; + double crtMatchingScore; }; // Internal data struct @@ -178,7 +183,8 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer { const geo::WireReadoutGeom *wireReadout, const detinfo::DetectorClocksData &dclock, const cheat::BackTrackerService *bt_serv, - const detinfo::DetectorPropertiesData &dprop); + const detinfo::DetectorPropertiesData &dprop, + const float &xshift); void FillTrackCRTHitInfo(const std::vector> &tag); // helpers @@ -192,6 +198,7 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer { art::InputTag fPFPT0producer; art::InputTag fCRTTrackT0producer; art::InputTag fCRTHitT0producer; + art::InputTag fCRTSpacePointT0producer; art::InputTag fCALOproducer; art::InputTag fTRKproducer; art::InputTag fTRKHMproducer; diff --git a/sbncode/Calibration/TrackCaloSkimmer_module.cc b/sbncode/Calibration/TrackCaloSkimmer_module.cc index 5fc33ab4d..6fe70bc1c 100644 --- a/sbncode/Calibration/TrackCaloSkimmer_module.cc +++ b/sbncode/Calibration/TrackCaloSkimmer_module.cc @@ -11,7 +11,6 @@ //////////////////////////////////////////////////////////////////////// #include "TrackCaloSkimmer.h" -#include "sbnobj/SBND/CRT/CRTTrack.hh" #include "art/Utilities/make_tool.h" @@ -70,6 +69,7 @@ sbn::TrackCaloSkimmer::TrackCaloSkimmer(fhicl::ParameterSet const& p) fPFPproducer = p.get< art::InputTag > ("PFPproducer","pandoraGausCryo0"); fPFPT0producer = p.get< art::InputTag > ("PFPT0producer", "pandoraGausCryo0"); fCRTTrackT0producer = p.get< art::InputTag >("CRTTrackT0producer", "crttrackmatching"); + fCRTSpacePointT0producer = p.get< art::InputTag >("CRTSpacePointT0producer", "crtspacepointmatching"); fCRTHitT0producer = p.get< art::InputTag >("CRTHitT0producer", "CRTT0Tagging"); fCALOproducer = p.get< art::InputTag > ("CALOproducer"); @@ -231,6 +231,8 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e) // // Tracks (SBND style) art::FindManyP fmT0CRTTrack(tracks, e, fCRTTrackT0producer); + // SpacePoints (SBND style) + art::FindManyP fmT0CRTSpacePoint(tracks, e, fCRTSpacePointT0producer); // Hits (ICARUS style) art::FindManyP fmT0CRTHit(tracks, e, fCRTHitT0producer); @@ -313,6 +315,9 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e) for (art::Ptr p_pfp: PFParticleList) { const recob::PFParticle &pfp = *p_pfp; + if(p_pfp->PdgCode() == 11) + continue; + const std::vector> thisTrack = fmTracks.at(p_pfp.key()); if (thisTrack.size() != 1) continue; @@ -345,8 +350,9 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e) } // Collect T0s - bool hasT0 = false, hasPFPT0 = false, hasCRTTrackT0 = false, hasCRTHitT0 = false; + bool hasT0 = false, hasPFPT0 = false, hasCRTTrackT0 = false, hasCRTHitT0 = false, hasCRTSpacePointT0 = false; int whicht0 = -1; + float crtMatchingScore = std::numeric_limits::signaling_NaN(); double t0PFP = std::numeric_limits::signaling_NaN(); if (fmT0PFP.isValid() && fmT0PFP.at(p_pfp.key()).size()) { @@ -356,12 +362,23 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e) } double t0CRTTrack = std::numeric_limits::signaling_NaN(); + double t0CRTTrackScore = std::numeric_limits::signaling_NaN(); if (fmT0CRTTrack.isValid() && fmT0CRTTrack.at(trkPtr.key()).size()) { t0CRTTrack = fmT0CRTTrack.data(trkPtr.key()).at(0)->Time(); + t0CRTTrackScore = fmT0CRTTrack.data(trkPtr.key()).at(0)->TriggerConfidence(); hasCRTTrackT0 = true; } + double t0CRTSpacePoint = std::numeric_limits::signaling_NaN(); + double t0CRTSpacePointScore = std::numeric_limits::signaling_NaN(); + if (fmT0CRTSpacePoint.isValid() && fmT0CRTSpacePoint.at(trkPtr.key()).size()) { + t0CRTSpacePoint = fmT0CRTSpacePoint.data(trkPtr.key()).at(0)->Time(); + t0CRTSpacePointScore = fmT0CRTSpacePoint.data(trkPtr.key()).at(0)->TriggerConfidence(); + hasCRTSpacePointT0 = true; + } + double t0CRTHit = std::numeric_limits::signaling_NaN(); + double t0CRTHitScore = std::numeric_limits::signaling_NaN(); if (fIncludeCRTHitTagging && fmT0CRTHit.isValid() && fmT0CRTHit.at(trkPtr.key()).size()) { const sbn::crt::CRTHitT0TaggingInfo &tag = *fmCRTHitT0TaggingInfo.at(trkPtr.key()).at(0); double time = fmT0CRTHit.at(trkPtr.key()).at(0)->Time(); @@ -384,23 +401,42 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e) if (!crtHitSysRejected && !crtHitDistanceRejected) { t0CRTHit = time; + t0CRTHitScore = tag.Distance; hasCRTHitT0 = true; } } - T0TimingInfo thisTrackTimingInfo = {t0PFP, t0CRTTrack, t0CRTHit, hasPFPT0, hasCRTTrackT0, hasCRTHitT0}; - hasT0 = hasPFPT0 || hasCRTTrackT0 || hasCRTHitT0; + hasT0 = hasPFPT0 || hasCRTTrackT0 || hasCRTHitT0 || hasCRTSpacePointT0; // "whicht0" should reflect the T0 used for the reconstruction of the drift coordinate. - if(!hasT0) whicht0 = -1 ; + if(!hasT0) whicht0 = -1; // In this way, if a track is T0 tagged from PFP and CRT tagged, which T0 reflects the PFP Tag. - else if (hasPFPT0) whicht0 = 0 ; - else if (hasCRTTrackT0) whicht0 = 1 ; - else if (hasCRTHitT0) whicht0 = 2 ; + else if (hasPFPT0) whicht0 = 0; + else if (hasCRTTrackT0) + { + whicht0 = 1; + crtMatchingScore = t0CRTTrackScore; + } + else if (hasCRTHitT0) + { + whicht0 = 2; + crtMatchingScore = t0CRTHitScore; + } + else if (hasCRTSpacePointT0) + { + whicht0 = 3; + crtMatchingScore = t0CRTSpacePointScore; + } + + T0TimingInfo thisTrackTimingInfo = {t0PFP, t0CRTTrack, t0CRTHit, t0CRTSpacePoint, hasPFPT0, hasCRTTrackT0, hasCRTHitT0, hasCRTSpacePointT0, crtMatchingScore}; if (fRequireT0 && !hasT0) continue; - if (fVerbose) std::cout << "Processing new track! ID: " << trkPtr->ID() << " time: " << t0PFP << " timeCRTTrack: " << t0CRTTrack << " timeCRTHit "<ID() << "\n" + << "\tPandora time: " << t0PFP << "\n" + << "\tCRTTrack (SBND) time: " << t0CRTTrack << "\n" + << "\tCRTHit (ICARUS) time: " << t0CRTHit << "\n" + << "\tCRTSpacePoint (SBND) time: " << t0CRTSpacePoint << std::endl; // Reset the track object *fTrack = sbn::TrackInfo(); @@ -1124,6 +1160,8 @@ void sbn::TrackCaloSkimmer::FillTrack(const recob::Track &track, fTrack->t0PFP = t0Info.t0Pandora; fTrack->t0CRTTrack = t0Info.t0CRTTrack; fTrack->t0CRTHit = t0Info.t0CRTHit; + fTrack->t0CRTSpacePoint = t0Info.t0CRTSpacePoint; + fTrack->crtMatchingScore = t0Info.crtMatchingScore; fTrack->id = track.ID(); fTrack->clear_cosmic_muon = pfp.Parent() == recob::PFParticle::kPFParticlePrimary; @@ -1135,22 +1173,40 @@ void sbn::TrackCaloSkimmer::FillTrack(const recob::Track &track, fTrack->dir.x = track.StartDirection().X(); fTrack->dir.y = track.StartDirection().Y(); fTrack->dir.z = track.StartDirection().Z(); - //If track is only CRTt0 tagged, undo the assumed trigger correction - if (t0Info.hasT0Pandora) { - fTrack->start.x = track.Start().X(); - fTrack->end.x = track.End().X(); - } else if (t0Info.hasT0CRTTrack) { - int driftDir = geo->TPC(hits[0]->WireID()).DriftDir().X(); - const double driftv(dprop.DriftVelocity(dprop.Efield(), dprop.Temperature())); - fTrack->start.x = track.Start().X() + driftDir*driftv*t0Info.t0CRTTrack*1e-3; - fTrack->end.x = track.End().X() + driftDir*driftv*t0Info.t0CRTTrack*1e-3; - } else if (t0Info.hasT0CRTHit){ - // If the track does not have a a Pandora T0, the tracks will always be either on the left or (ex Or) right of the cathode. - int driftDir = geo->TPC(hits[0]->WireID()).DriftDir().X(); - const double driftv(dprop.DriftVelocity(dprop.Efield(), dprop.Temperature())); - fTrack->start.x = track.Start().X() + driftDir*driftv*t0Info.t0CRTHit*1e-3; - fTrack->end.x = track.End().X() + driftDir*driftv*t0Info.t0CRTHit*1e-3; - } + + if(t0Info.hasT0Pandora) + { + fTrack->start.x = track.Start().X(); + fTrack->end.x = track.End().X(); + } + else + { + // If CRT T0 tagged we can shift the x positions appropriately + int driftDir = geo->TPC(hits[0]->WireID()).DriftDir().X(); + const double driftv(dprop.DriftVelocity(dprop.Efield(), dprop.Temperature())); + + if(t0Info.hasT0CRTTrack) + { + const double xshift = driftDir*driftv*t0Info.t0CRTTrack*1e-3; + fTrack->start.x = track.Start().X() + xshift; + fTrack->end.x = track.End().X() + xshift; + fTrack->xShiftCRT = xshift; + } + else if(t0Info.hasT0CRTHit) + { + const double xshift = driftDir*driftv*t0Info.t0CRTHit*1e-3; + fTrack->start.x = track.Start().X() + xshift; + fTrack->end.x = track.End().X() + xshift; + fTrack->xShiftCRT = xshift; + } + else if(t0Info.hasT0CRTSpacePoint) + { + const double xshift = driftDir*driftv*t0Info.t0CRTSpacePoint*1e-3; + fTrack->start.x = track.Start().X() + xshift; + fTrack->end.x = track.End().X() + xshift; + fTrack->xShiftCRT = xshift; + } + } if (hits.size() > 0) { fTrack->cryostat = hits[0]->WireID().Cryostat; @@ -1158,7 +1214,7 @@ void sbn::TrackCaloSkimmer::FillTrack(const recob::Track &track, // Fill each hit for (unsigned i_hit = 0; i_hit < hits.size(); i_hit++) { - sbn::TrackHitInfo hinfo = MakeHit(*hits[i_hit], hits[i_hit].key(), *thms[i_hit], track, t0Info, sps[i_hit], calo, geo, wireReadout, clock_data, bt_serv, dprop); + sbn::TrackHitInfo hinfo = MakeHit(*hits[i_hit], hits[i_hit].key(), *thms[i_hit], track, t0Info, sps[i_hit], calo, geo, wireReadout, clock_data, bt_serv, dprop, fTrack->xShiftCRT); if (hinfo.h.plane == 0) { fTrack->hits0.push_back(hinfo); } @@ -1325,7 +1381,8 @@ sbn::TrackHitInfo sbn::TrackCaloSkimmer::MakeHit(const recob::Hit &hit, const geo::WireReadoutGeom *wireReadout, const detinfo::DetectorClocksData &dclock, const cheat::BackTrackerService *bt_serv, - const detinfo::DetectorPropertiesData &dprop) { + const detinfo::DetectorPropertiesData &dprop, + const float &xshift) { // TrackHitInfo to save sbn::TrackHitInfo hinfo; @@ -1400,13 +1457,29 @@ sbn::TrackHitInfo sbn::TrackCaloSkimmer::MakeHit(const recob::Hit &hit, } } - // This is needed to reconstrut drift coordinate using a different Time - int driftDir = geo->TPC(hit.WireID()).DriftDir().X(); - const double driftv(dprop.DriftVelocity(dprop.Efield(), dprop.Temperature())); - double anodeDistance = std::numeric_limits::signaling_NaN(); - if(t0Info.hasT0CRTHit) anodeDistance = (hit.PeakTime()-dclock.Time2Tick(dclock.TriggerTime())-t0Info.t0CRTHit*1e-3/dclock.TPCClock().TickPeriod())*dclock.TPCClock().TickPeriod()*driftv; - double wirePlaneX = wireReadout->Plane(geo::PlaneID(hit.WireID().Cryostat, hit.WireID().TPC, hit.WireID().Plane)).GetCenter().X(); - double recoX = wirePlaneX - driftDir*anodeDistance; + // If T0 provided CRT then we should appropriately shift the drift position (already done for pandora T0) + double recoX = std::numeric_limits::signaling_NaN(); + + if(!t0Info.hasT0Pandora && (t0Info.hasT0CRTTrack || t0Info.hasT0CRTHit || t0Info.hasT0CRTSpacePoint)) + { + int driftDir = geo->TPC(hit.WireID()).DriftDir().X(); + const double driftv(dprop.DriftVelocity(dprop.Efield(), dprop.Temperature())); + + double time = std::numeric_limits::signaling_NaN(); + + if(t0Info.hasT0CRTTrack) + time = t0Info.t0CRTTrack*1e-3; + else if(t0Info.hasT0CRTHit) + time = t0Info.t0CRTHit*1e-3; + else if(t0Info.hasT0CRTSpacePoint) + time = t0Info.t0CRTSpacePoint*1e-3; + + double anodeDistance = (hit.PeakTime()-dclock.Time2Tick(dclock.TriggerTime())-time/dclock.TPCClock().TickPeriod())*dclock.TPCClock().TickPeriod()*driftv; + double wirePlaneX = wireReadout->Plane(geo::PlaneID(hit.WireID().Cryostat, hit.WireID().TPC, hit.WireID().Plane)).GetCenter().X(); + + recoX = wirePlaneX - driftDir*anodeDistance; + } + // Information from the TrackHitMeta bool badhit = (thm.Index() == std::numeric_limits::max()) || (!trk.HasValidPoint(thm.Index())); @@ -1419,7 +1492,7 @@ sbn::TrackHitInfo sbn::TrackCaloSkimmer::MakeHit(const recob::Hit &hit, // The tp X coordinate is reconstructed only if it does not have a PandoraT0. hinfo.tp.x = loc.X(); - if(t0Info.hasT0CRTHit && !t0Info.hasT0Pandora && !isnan(hinfo.tp.x)) hinfo.tp.x = recoX; + if((t0Info.hasT0CRTTrack || t0Info.hasT0CRTHit || t0Info.hasT0CRTSpacePoint) && !t0Info.hasT0Pandora && !isnan(hinfo.tp.x)) hinfo.tp.x = recoX; hinfo.tp.y = loc.Y(); hinfo.tp.z = loc.Z(); @@ -1451,7 +1524,7 @@ sbn::TrackHitInfo sbn::TrackCaloSkimmer::MakeHit(const recob::Hit &hit, if (sp) { hinfo.h.sp.x = sp->position().x(); // If the track has a Pandora T0, do not displace, track is already reconstructed at the correct position - if(t0Info.hasT0CRTHit && !t0Info.hasT0Pandora) hinfo.h.sp.x = sp->position().x() + driftDir*driftv*t0Info.t0CRTHit*1e-3; + if((t0Info.hasT0CRTTrack || t0Info.hasT0CRTHit || t0Info.hasT0CRTSpacePoint) && !t0Info.hasT0Pandora) hinfo.h.sp.x = sp->position().x() + xshift; hinfo.h.sp.y = sp->position().y(); hinfo.h.sp.z = sp->position().z(); hinfo.h.hasSP = true; diff --git a/sbncode/Calibration/fcl/sbnd/sbnd_trackcalo_skimmer.fcl b/sbncode/Calibration/fcl/sbnd/sbnd_trackcalo_skimmer.fcl index 8735ed405..b510a231f 100644 --- a/sbncode/Calibration/fcl/sbnd/sbnd_trackcalo_skimmer.fcl +++ b/sbncode/Calibration/fcl/sbnd/sbnd_trackcalo_skimmer.fcl @@ -83,6 +83,6 @@ a2c_selection: { } caloskim_nodigits_goldentracks: @local::caloskim_nodigits -caloskim_nodigits_goldentracks.SelectionTools: [@local::stopping_selection, @local::a2c_selection, @local::throughgoing_selection] +caloskim_nodigits_goldentracks.SelectionTools: [ @local::stopping_selection, @local::a2c_selection, @local::throughgoing_selection ] END_PROLOG