26 #include "canvas/Persistency/Common/FindManyP.h" 27 #include "art_root_io/TFileService.h" 35 #include "nug4/ParticleNavigation/ParticleList.h" 74 class SingleCRTMatchingProducer;
91 bool moduleMatcher(
int module1,
int module2);
113 double trackX1, trackX2, trackY1, trackY2, trackZ1,
trackZ2;
201 consumes < std::vector < CRT::Trigger >> (
fCRTLabel);
202 consumes < std::vector < art::Assns < sim::AuxDetSimChannel, CRT::Trigger >>> (
fCRTLabel);
204 produces< std::vector<anab::T0> >();
205 produces< std::vector<anab::CosmicTag> >();
206 produces< art::Assns<recob::Track, anab::T0> >();
207 produces< art::Assns<recob::Track, anab::CosmicTag> >();
208 produces< art::Assns<anab::CosmicTag, anab::T0> >();
209 produces< art::Assns<CRT::Trigger, anab::CosmicTag> >();
218 if ((module1 == 6 && (module2 == 10 || module2 == 11)) || (module1 == 14 && (module2 == 10 || module2 == 11)) || (module1 == 19 && (module2 == 26 || module2 == 27)) || (module1 == 31 && (module2 == 26 || module2 == 27)) || (module1 == 7 && (module2 == 12 || module2 == 13)) || (module1 == 15 && (module2 == 12 || module2 == 13)) || (module1 == 18 && (module2 == 24 || module2 == 25)) || (module1 == 30 && (module2 == 24 || module2 == 25)) || (module1 == 1 && (module2 == 4 || module2 == 5)) || (module1 == 9 && (module2 == 4 || module2 == 5)) || (module1 == 16 && (module2 == 20 || module2 == 21)) || (module1 == 28 && (module2 == 20 || module2 == 21)) || (module1 == 0 && (module2 == 2 || module2 == 3)) || (module1 == 8 && (module2 == 2 || module2 == 3)) || (module1 == 17 && (module2 == 22 || module2 == 23)) || (module1 == 29 && (module2 == 22 || module2 == 23)))
return 1;
228 std::unique_ptr< std::vector<anab::T0> > T0col(
new std::vector<anab::T0>);
230 auto CRTTrack=std::make_unique< std::vector< anab::CosmicTag > > ();
281 const auto & triggers =
event.getValidHandle < std::vector < CRT::Trigger >> (
fCRTLabel);
283 art::FindManyP < sim::AuxDetSimChannel > trigToSim(triggers, event,
fCRTLabel);
285 vector < art::Ptr < CRT::Trigger > > crtList;
286 auto crtListHandle =
event.getHandle < vector < CRT::Trigger > >(
fCRTLabel);
293 auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
296 std::unordered_map < size_t, double > prevTimes;
299 for (
const auto &
trigger: * triggers) {
300 const auto & hits =
trigger.Hits();
301 for (
const auto &
hit: hits) {
326 const auto & trigGeo = geom -> AuxDet(
trigger.Channel());
327 const auto & csens = trigGeo.SensitiveVolume(
hit.Channel());
328 const auto center = csens.GetCenter();
341 for (
unsigned int f_test = 0; f_test <
tempHits_F.size(); f_test++) {
343 const auto & trigGeo = geom -> AuxDet(
tempHits_F[
f].module);
344 const auto & trigGeo2 = geom -> AuxDet(
tempHits_F[f_test].module);
346 const auto & hit1Geo = trigGeo.SensitiveVolume(
tempHits_F[
f].channelGeo);
347 const auto hit1Center = hit1Geo.GetCenter();
349 const auto & hit2Geo = trigGeo2.SensitiveVolume(
tempHits_F[f_test].channelGeo);
350 const auto hit2Center = hit2Geo.GetCenter();
356 double hitX = hit1Center.X();
361 double hitYPrelim=hit2Center.Y();
369 double hitY=hitYPrelim;
370 double hitZ = (hit1Center.Z() + hit2Center.Z()) / 2.
f;
390 for (
unsigned int f_test = 0; f_test <
tempHits_B.size(); f_test++) {
393 const auto & trigGeo = geom -> AuxDet(
tempHits_B[
f].module);
394 const auto & trigGeo2 = geom -> AuxDet(
tempHits_B[f_test].module);
395 const auto & hit1Geo = trigGeo.SensitiveVolume(
tempHits_B[
f].channelGeo);
396 const auto hit1Center = hit1Geo.GetCenter();
398 const auto & hit2Geo = trigGeo2.SensitiveVolume(
tempHits_B[f_test].channelGeo);
399 const auto hit2Center = hit2Geo.GetCenter();
405 double hitX = hit1Center.X();
413 double hitYPrelim = hit2Center.Y();
419 double hitY=hitYPrelim;
422 double hitZ = (hit1Center.Z() + hit2Center.Z()) / 2.
f;
448 vector < art::Ptr < recob::Track > > trackList;
449 auto trackListHandle =
event.getHandle < vector < recob::Track > >(
fTrackModuleLabel);
450 if (trackListHandle) {
454 vector<art::Ptr<recob::PFParticle> > pfplist;
455 auto PFPListHandle =
event.getHandle< std::vector<recob::PFParticle> >(
"pandora");
457 if (pfplist.size()<1)
return;
458 art::FindManyP<anab::T0> trk_t0_assn_v(PFPListHandle, event ,
"pandora");
459 art::FindManyP<recob::PFParticle> pfp_trk_assn(trackListHandle,event,
"pandoraTrack");
461 int nTracksReco = trackList.size();
462 art::FindManyP < recob::Hit > hitsFromTrack(trackListHandle, event,
fTrackModuleLabel);
467 for (
int iRecoTrack = 0; iRecoTrack < nTracksReco; ++iRecoTrack) {
469 std::vector< art::Ptr<recob::Hit> > allHits = hitsFromTrack.at(iRecoTrack);
473 std::vector<art::Ptr<recob::PFParticle>> pfps=pfp_trk_assn.at(iRecoTrack);
474 if(!pfps.size())
continue;
475 std::vector<art::Ptr<anab::T0>> t0s=trk_t0_assn_v.at(pfps[0].
key());
482 int lastHit=allHits.size()-2;
487 double trackStartPositionZ_noSCE = trackList[iRecoTrack]->Vertex().Z();
488 double trackEndPositionZ_noSCE = trackList[iRecoTrack] ->
End().Z();
490 double trackStartPositionX_noSCE = trackList[iRecoTrack]->Vertex().X();
491 double trackStartPositionY_noSCE = trackList[iRecoTrack]->Vertex().Y();
494 double trackEndPositionX_noSCE = trackList[iRecoTrack] ->
End().X();
495 double trackEndPositionY_noSCE = trackList[iRecoTrack] ->
End().Y();
497 if (trackStartPositionZ_noSCE>trackEndPositionZ_noSCE){
498 trackEndPositionZ_noSCE = trackList[iRecoTrack]->Vertex().Z();
499 trackStartPositionZ_noSCE = trackList[iRecoTrack] ->
End().Z();
500 trackEndPositionX_noSCE = trackList[iRecoTrack]->Vertex().X();
501 trackEndPositionY_noSCE = trackList[iRecoTrack]->Vertex().Y();
504 trackStartPositionX_noSCE = trackList[iRecoTrack] ->
End().X();
505 trackStartPositionY_noSCE = trackList[iRecoTrack] ->
End().Y();
509 if ((trackEndPositionZ_noSCE>90 && trackEndPositionZ_noSCE < 660 && trackStartPositionZ_noSCE <50 && trackStartPositionZ_noSCE<660) || (trackStartPositionZ_noSCE>90 && trackStartPositionZ_noSCE < 660 && trackEndPositionZ_noSCE <50 && trackEndPositionZ_noSCE<660)) {
511 double min_delta = DBL_MAX;
513 double best_dotProductCos = DBL_MAX;
514 double best_deltaXF = DBL_MAX;
515 double best_deltaYF = DBL_MAX;
516 int bestHitIndex_F=0;
517 double best_trackX1=DBL_MAX;
518 double best_trackX2=DBL_MAX;
519 for (
unsigned int iHit_F = 0; iHit_F <
primaryHits_F.size(); iHit_F++) {
522 double trackStartPositionX_notCorrected=trackStartPositionX_noSCE;
523 double trackEndPositionX_notCorrected=trackEndPositionX_noSCE;
531 double ticksOffset=0;
534 if (!
fMCCSwitch) ticksOffset = (
primaryHits_F[iHit_F].timeAvg+RDOffset)/25.
f+detProp.GetXTicksOffset(allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
536 else if (
fMCCSwitch) ticksOffset = (
primaryHits_F[iHit_F].timeAvg/500.f)+detProp.GetXTicksOffset(allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
538 xOffset=detProp.ConvertTicksToX(ticksOffset,allHits[firstHit]->
WireID().
Plane, allHits[firstHit]->
WireID().
TPC, allHits[firstHit]->
WireID().Cryostat);
541 trackStartPositionX_noSCE=trackStartPositionX_notCorrected-xOffset;
542 trackEndPositionX_noSCE=trackEndPositionX_notCorrected-xOffset;
543 if (fabs(xOffset)>300 || ((trackStartPositionX_notCorrected<0 && trackStartPositionX_noSCE>0) || (trackEndPositionX_notCorrected<0 && trackEndPositionX_noSCE>0)) || ((trackStartPositionX_notCorrected>0 && trackStartPositionX_noSCE<0) || (trackEndPositionX_notCorrected>0 && trackEndPositionX_noSCE<0)) )
continue;
546 double trackStartPositionX=trackStartPositionX_noSCE;
547 double trackStartPositionY=trackStartPositionY_noSCE;
548 double trackStartPositionZ=trackStartPositionZ_noSCE;
550 double trackEndPositionX=trackEndPositionX_noSCE;
551 double trackEndPositionY=trackEndPositionY_noSCE;
552 double trackEndPositionZ=trackEndPositionZ_noSCE;
557 trackStartPositionX -= posOffsets_F.X();
558 trackStartPositionY += posOffsets_F.Y();
559 trackStartPositionZ += posOffsets_F.Z();
561 trackEndPositionX -= posOffsets_B.X();
562 trackEndPositionY += posOffsets_B.Y();
563 trackEndPositionZ += posOffsets_B.Z();
574 TVector3 trackStart(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
575 TVector3 trackEnd(trackEndPositionX, trackEndPositionY, trackEndPositionZ);
576 TVector3 v1(X1,Y1,Z1);
577 TVector3 v2(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
579 TVector3 v4(trackStartPositionX,
581 trackStartPositionZ);
582 TVector3 v5(trackEndPositionX,
585 TVector3 trackVector = (v5-v4).Unit();
586 TVector3 hitVector=(v2-v1).Unit();
591 double predictedHitPositionY1 = (v1.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.Y()-v5.Y())+v5.Y();
594 double predictedHitPositionX1 = (v1.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.X()-v5.X())+v5.X();
597 if (predictedHitPositionX1<-200 || predictedHitPositionX1>580 || predictedHitPositionY1<-50 || predictedHitPositionY1>620)
continue;
599 double dotProductCos=trackVector*hitVector;
601 double deltaX1 = (predictedHitPositionX1-X1);
604 double deltaY1 = (predictedHitPositionY1-Y1);
611 best_dotProductCos = dotProductCos;
612 best_trackX1=trackStartPositionX;
613 best_trackX2=trackEndPositionX;
614 best_deltaXF = deltaX1;
615 best_deltaYF = deltaY1;
616 bestHitIndex_F=iHit_F;
620 int f=bestHitIndex_F;
623 double dotProductCos=best_dotProductCos;
625 double trackStartPositionX=best_trackX1;
626 double trackStartPositionY=trackStartPositionY_noSCE;
627 double trackStartPositionZ=trackStartPositionZ_noSCE;
629 double trackEndPositionX=best_trackX2;
630 double trackEndPositionY=trackEndPositionY_noSCE;
631 double trackEndPositionZ=trackEndPositionZ_noSCE;
634 TVector3 trackStart(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
635 TVector3 trackEnd(trackEndPositionX, trackEndPositionY, trackEndPositionZ);
640 tPair.
recoId = iRecoTrack;
666 if ( (trackStartPositionZ_noSCE<620 && trackEndPositionZ_noSCE > 660 && trackStartPositionZ_noSCE > 50 && trackEndPositionZ_noSCE > 50) || (trackStartPositionZ_noSCE>660 && trackEndPositionZ_noSCE < 620 && trackStartPositionZ_noSCE > 50 && trackEndPositionZ_noSCE > 50)) {
667 double min_delta = DBL_MAX;
669 double best_dotProductCos = DBL_MAX;
670 double best_deltaXF = DBL_MAX;
671 double best_deltaYF = DBL_MAX;
672 int bestHitIndex_B=0;
673 double best_trackX1=DBL_MAX;
674 double best_trackX2=DBL_MAX;
675 for (
unsigned int iHit_B = 0; iHit_B <
primaryHits_B.size(); iHit_B++) {
679 double trackStartPositionX_notCorrected=trackStartPositionX_noSCE;
680 double trackEndPositionX_notCorrected=trackEndPositionX_noSCE;
688 double ticksOffset=0;
689 if (!
fMCCSwitch) ticksOffset = (
primaryHits_B[iHit_B].timeAvg+RDOffset)/25.
f+detProp.GetXTicksOffset(allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
691 else if (
fMCCSwitch) ticksOffset = (
primaryHits_B[iHit_B].timeAvg/500.f)+detProp.GetXTicksOffset(allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
693 xOffset=detProp.ConvertTicksToX(ticksOffset,allHits[firstHit]->
WireID().
Plane, allHits[firstHit]->
WireID().
TPC, allHits[firstHit]->
WireID().Cryostat);
696 trackStartPositionX_noSCE=trackStartPositionX_notCorrected-xOffset;
697 trackEndPositionX_noSCE=trackEndPositionX_notCorrected-xOffset;
699 if (fabs(xOffset)>300 || ((trackStartPositionX_notCorrected<0 && trackStartPositionX_noSCE>0) || (trackEndPositionX_notCorrected<0 && trackEndPositionX_noSCE>0)) || ((trackStartPositionX_notCorrected>0 && trackStartPositionX_noSCE<0) || (trackEndPositionX_notCorrected>0 && trackEndPositionX_noSCE<0)) )
continue;
702 double trackStartPositionX=trackStartPositionX_noSCE;
703 double trackStartPositionY=trackStartPositionY_noSCE;
704 double trackStartPositionZ=trackStartPositionZ_noSCE;
706 double trackEndPositionX=trackEndPositionX_noSCE;
707 double trackEndPositionY=trackEndPositionY_noSCE;
708 double trackEndPositionZ=trackEndPositionZ_noSCE;
718 TVector3 trackStart(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
719 TVector3 trackEnd(trackEndPositionX, trackEndPositionY, trackEndPositionZ);
720 TVector3 v1(X1,Y1,Z1);
721 TVector3 v2(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
723 TVector3 v4(trackStartPositionX,
725 trackStartPositionZ);
726 TVector3 v5(trackEndPositionX,
729 TVector3 trackVector = (v5-v4).Unit();
730 TVector3 hitVector=(v2-v1).Unit();
735 double predictedHitPositionY1 = (v1.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.Y()-v5.Y())+v5.Y();
738 double predictedHitPositionX1 = (v1.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.X()-v5.X())+v5.X();
740 if (
abs(predictedHitPositionX1)>340 || predictedHitPositionY1<-160 || predictedHitPositionY1>560)
continue;
742 double dotProductCos=trackVector*hitVector;
744 double deltaX1 = (predictedHitPositionX1-X1);
746 double deltaY1 = (predictedHitPositionY1-Y1);
752 best_dotProductCos = dotProductCos;
753 best_trackX1=trackStartPositionX;
754 best_trackX2=trackEndPositionX;
755 best_deltaXF = deltaX1;
756 best_deltaYF = deltaY1;
757 bestHitIndex_B=iHit_B;
761 int f=bestHitIndex_B;
764 double dotProductCos=best_dotProductCos;
767 double trackStartPositionX=best_trackX1;
768 double trackStartPositionY=trackStartPositionY_noSCE;
769 double trackStartPositionZ=trackStartPositionZ_noSCE;
771 double trackEndPositionX=best_trackX2;
772 double trackEndPositionY=trackEndPositionY_noSCE;
773 double trackEndPositionZ=trackEndPositionZ_noSCE;
775 TVector3 trackStart(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
776 TVector3 trackEnd(trackEndPositionX, trackEndPositionY, trackEndPositionZ);
782 tPair.
recoId = iRecoTrack;
816 vector < tracksPair > allUniqueTracksPair;
830 if (allUniqueTracksPair.size() > 0) {
831 for (
unsigned int u = 0; u < allUniqueTracksPair.size(); u++) {
833 int CRTTrackId=allUniqueTracksPair[u].CRTTrackId;
834 int TPCTrackId=allUniqueTracksPair[u].recoId;
835 deltaX=allUniqueTracksPair[u].deltaX;
837 deltaY=allUniqueTracksPair[u].deltaY;
840 dotCos=fabs(allUniqueTracksPair[u].dotProductCos);
841 trackX1=allUniqueTracksPair[u].trackStartPosition.X();
842 trackY1=allUniqueTracksPair[u].trackStartPosition.Y();
843 trackZ1=allUniqueTracksPair[u].trackStartPosition.Z();
845 trackX2=allUniqueTracksPair[u].trackEndPosition.X();
846 trackY2=allUniqueTracksPair[u].trackEndPosition.Y();
847 trackZ2=allUniqueTracksPair[u].trackEndPosition.Z();
849 moduleX=allUniqueTracksPair[u].moduleX1;
850 moduleY=allUniqueTracksPair[u].moduleY1;
852 adcX=allUniqueTracksPair[u].adcX1;
853 adcY=allUniqueTracksPair[u].adcY1;
855 CRTT0=allUniqueTracksPair[u].timeAvg;
859 stripX=allUniqueTracksPair[u].stripX1;
860 stripY=allUniqueTracksPair[u].stripY1;
862 X_CRT=allUniqueTracksPair[u].X1;
863 Y_CRT=allUniqueTracksPair[u].Y1;
864 Z_CRT=allUniqueTracksPair[u].Z1;
871 std::vector<float> hitF;
872 std::vector<float> hitB;
876 if (
Z_CRT<100) T0col->push_back(
anab::T0(CRTT0, 1, 1,CRTTrackId,fabs(allUniqueTracksPair[u].dotProductCos)));
877 else T0col->push_back(
anab::T0(CRTT0, 2, 1,CRTTrackId,fabs(allUniqueTracksPair[u].dotProductCos) ));
878 auto const crtTrackPtr = crtPtr(CRTTrack->size()-1);
879 auto const t0CP = t0CandPtr(CRTTrackId);
880 CRTT0assn->addSingle(crtTrackPtr,t0CP);
883 util::CreateAssn(*
this, event, *CRTTrack, trackList[TPCTrackId], *TPCCRTassn);
884 util::CreateAssn(*
this, event, *CRTTrack, crtList[allUniqueTracksPair[u].trigNumberX], *CRTTriggerassn);
885 util::CreateAssn(*
this, event, *CRTTrack, crtList[allUniqueTracksPair[u].trigNumberY], *CRTTriggerassn);
int fModuletoModuleTimingCut
std::vector< tempHits > tempHits_F
std::string fTrackModuleLabel
TVector3 trackStartPosition
geo::TPCID PositionToTPCID(geo::Point_t const &point) const
Returns the ID of the TPC at specified location.
std::vector< recoHits > primaryHits_F
std::vector< recoHits > primaryHits_B
SingleCRTMatchingProducer(fhicl::ParameterSet const &p)
art framework interface to geometry description
std::vector< tracksPair > tracksPair_B
TVector3 trackEndPosition
#define DEFINE_ART_MODULE(klass)
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
void produce(art::Event &event) override
removePairIndex(const tracksPair &tracksPair0)
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
bool operator()(const tracksPair &tracksPair2)
std::vector< tempHits > tempHits_B
Detector simulation of raw signals on wires.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
const tracksPair tracksPair1
Declaration of signal hit object.
bool operator()(const tracksPair &pair1, const tracksPair &pair2)
int fFronttoBackTimingCut
Provides recob::Track data product.
bool moduleMatcher(int module1, int module2)
constexpr auto const & deepestIndex() const
Returns the value of the deepest ID available (TPC's).
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
recob::tracking::Plane Plane
std::vector< tracksPair > tracksPair_F
Event finding and building.