Inheritance diagram for CRT::TwoCRTMatchingProducer:

Classes
struct	recoHits

struct	tempHits

Public Member Functions
	TwoCRTMatchingProducer (fhicl::ParameterSet const &p)

	TwoCRTMatchingProducer (TwoCRTMatchingProducer const &)=delete

	TwoCRTMatchingProducer (TwoCRTMatchingProducer &&)=delete

TwoCRTMatchingProducer &	operator= (TwoCRTMatchingProducer const &)=delete

TwoCRTMatchingProducer &	operator= (TwoCRTMatchingProducer &&)=delete

bool	moduleMatcher (int module1, int module2)

void	produce (art::Event &event) override

Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)

template<typename Config >
	EDProducer (Table< Config > const &config)

std::string	workerType () const

Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept

	Producer (fhicl::ParameterSet const &)

	Producer (Producer const &)=delete

	Producer (Producer &&)=delete

Producer &	operator= (Producer const &)=delete

Producer &	operator= (Producer &&)=delete

void	doBeginJob (SharedResources const &resources)

void	doEndJob ()

void	doRespondToOpenInputFile (FileBlock const &fb)

void	doRespondToCloseInputFile (FileBlock const &fb)

void	doRespondToOpenOutputFiles (FileBlock const &fb)

void	doRespondToCloseOutputFiles (FileBlock const &fb)

bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)

Public Member Functions inherited from art::Modifier
	~Modifier () noexcept

	Modifier ()

	Modifier (Modifier const &)=delete

	Modifier (Modifier &&)=delete

Modifier &	operator= (Modifier const &)=delete

Modifier &	operator= (Modifier &&)=delete

Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept

	ModuleBase ()

ModuleDescription const &	moduleDescription () const

void	setModuleDescription (ModuleDescription const &)

std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const

void	sortConsumables (std::string const &current_process_name)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Public Attributes
int	nEvents = 0

int	nHitsPerEvent =0

std::string	fTrackModuleLabel = "pandoraTrack"

Private Attributes
art::InputTag	fCRTLabel

bool	fSCECorrection

bool	fModuleSwitch

int	fADCThreshold

int	fModuletoModuleTimingCut

int	fFronttoBackTimingCut

long long	timeStamp

std::vector< recoHits >	primaryHits_F

std::vector< recoHits >	primaryHits_B

std::vector< tempHits >	tempHits_F

std::vector< tempHits >	tempHits_B

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >

Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >

Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)

Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 85 of file TwoCRTMatchingProducer_module.cc.

Constructor & Destructor Documentation

CRT::TwoCRTMatchingProducer::TwoCRTMatchingProducer ( fhicl::ParameterSet const & p )

explicit

Definition at line 149 of file TwoCRTMatchingProducer_module.cc.

                                                                             : EDProducer{p}, fCRTLabel(p.get < art::InputTag > ("CRTLabel"))
 {
   consumes < std::vector < CRT::Trigger >> (fCRTLabel);
   consumes < std::vector < art::Assns < sim::AuxDetSimChannel, CRT::Trigger >>> (fCRTLabel);
   produces< std::vector<anab::T0> >();
 
   produces< std::vector<anab::CosmicTag> >();  
   produces< art::Assns<recob::Track, anab::T0> >();
   produces< art::Assns<recob::Track, anab::CosmicTag> >();
   produces< art::Assns<CRT::Trigger, anab::CosmicTag> >();
 
   fSCECorrection=(p.get<bool>("SCECorrection"));
 }

CRT::TwoCRTMatchingProducer::TwoCRTMatchingProducer ( TwoCRTMatchingProducer const & )

delete

CRT::TwoCRTMatchingProducer::TwoCRTMatchingProducer ( TwoCRTMatchingProducer && )

delete

Member Function Documentation

bool CRT::TwoCRTMatchingProducer::moduleMatcher	(	int	module1,
		int	module2
	)

Definition at line 164 of file TwoCRTMatchingProducer_module.cc.

                                                                       {
   // Function checking if two hits could reasonably be matched into a 2D hit
   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;
   else return 0;
 
 }

TwoCRTMatchingProducer& CRT::TwoCRTMatchingProducer::operator= ( TwoCRTMatchingProducer const & )

delete

TwoCRTMatchingProducer& CRT::TwoCRTMatchingProducer::operator= ( TwoCRTMatchingProducer && )

delete

void CRT::TwoCRTMatchingProducer::produce ( art::Event & event )

overridevirtual

Implements art::EDProducer.

Definition at line 173 of file TwoCRTMatchingProducer_module.cc.

 {
     bool fMCCSwitch=!event.isRealData();
     nEvents++;  
   std::unique_ptr< std::vector<anab::T0> > T0col( new std::vector<anab::T0>);
   auto CRTTrack=std::make_unique< std::vector< anab::CosmicTag > > (); 
 
 
  std::unique_ptr< art::Assns<CRT::Trigger, anab::CosmicTag>> CRTTriggerassn( new art::Assns<CRT::Trigger, anab::CosmicTag>);
 
  std::unique_ptr< art::Assns<recob::Track, anab::CosmicTag> > TPCCRTassn( new art::Assns<recob::Track, anab::CosmicTag>);
  std::unique_ptr< art::Assns<recob::Track, anab::T0> > TPCT0assn( new art::Assns<recob::Track, anab::T0>);
    
 
   if (fMCCSwitch){
     fModuleSwitch=1;
     fADCThreshold=800;
     fModuletoModuleTimingCut=2;
     fFronttoBackTimingCut=100;
     
 }
   else {
     fModuleSwitch=0;
     fADCThreshold=10;
     fModuletoModuleTimingCut=5;
     fFronttoBackTimingCut=8;
     art::ValidHandle<std::vector<raw::RDTimeStamp>> timingHandle = event.getValidHandle<std::vector<raw::RDTimeStamp>>("timingrawdecoder:daq");
     timeStamp=timingHandle->at(0).GetTimeStamp();
 } 
   int nHits = 0;
 
         //Detector properties service
   auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(event);
 
   primaryHits_F.clear();
   primaryHits_B.clear();
   tempHits_F.clear();
   tempHits_B.clear(); // Arrays to compile hits and move them through
 
 
   //Get triggers
   //cout << "Getting triggers" << endl;
   vector < art::Ptr < CRT::Trigger > > crtList;
   auto crtListHandle = event.getHandle < vector < CRT::Trigger > >(fCRTLabel);
   if (crtListHandle) {
     art::fill_ptr_vector(crtList, crtListHandle);
   }
   const auto & triggers = event.getValidHandle < std::vector < CRT::Trigger >> (fCRTLabel);
 
   art::FindManyP < sim::AuxDetSimChannel > trigToSim(triggers, event, fCRTLabel);
 
   //Get a handle to the Geometry service to look up AuxDetGeos from module numbers
   art::ServiceHandle < geo::Geometry > geom;
 
   auto const* SCE = lar::providerFrom<spacecharge::SpaceChargeService>();
 
   //Mapping from channel to trigger
   std::unordered_map < size_t, double > prevTimes;
   int hitID = 0;
   //cout << "Looking for hits in Triggers" << endl;
   int trigID=0;
   for (const auto & trigger: * triggers) {
     const auto & hits = trigger.Hits();
     for (const auto & hit: hits) { // Collect hits on all modules
         //cout<<hits.size()<<','<<hit.ADC()<<endl;
       if (hit.ADC() > fADCThreshold) { // Keep if they are above threshold
 
         tempHits tHits;
         if (!fMCCSwitch){
 
         tHits.module = trigger.Channel(); // Values to add to array
         tHits.channel=hit.Channel();
         tHits.adc = hit.ADC();
         tHits.triggerTime=trigger.Timestamp()-timeStamp;
         }
         else{
         tHits.module = trigger.Channel(); // Values to add to array
         tHits.channel=hit.Channel();
         tHits.adc = hit.ADC();
         tHits.triggerTime=trigger.Timestamp();
         }
          //cout<<trigger.Channel()<<','<<hit.Channel()<<','<<hit.ADC()<<endl;
         nHits++;
         tHits.triggerNumber=trigID;
         const auto & trigGeo = geom -> AuxDet(trigger.Channel()); // Get geo  
         const auto & csens = trigGeo.SensitiveVolume(hit.Channel());
         const auto center = csens.GetCenter();
         if (center.Z() < 100) tempHits_F.push_back(tHits); // Sort F/B from Z
         else tempHits_B.push_back(tHits);
         hitID++;
       }
     }
     trigID++;
   }
   nHitsPerEvent=nHits;
   //cout << "Hits compiled for event: " << nEvents << endl;
   //cout << "Number of Hits above Threshold:  " << hitID << endl;
 
   for (unsigned int f = 0; f < tempHits_F.size(); f++) {
     for (unsigned int f_test = 0; f_test < tempHits_F.size(); f_test++) {
        if (fabs(tempHits_F[f_test].triggerTime-tempHits_F[f].triggerTime)>fModuletoModuleTimingCut) continue;
       const auto & trigGeo = geom -> AuxDet(tempHits_F[f].module);
       const auto & trigGeo2 = geom -> AuxDet(tempHits_F[f_test].module);
 
       const auto & hit1Geo = trigGeo.SensitiveVolume(tempHits_F[f].channel);
       const auto hit1Center = hit1Geo.GetCenter();
       // Create 2D hits from geo of the Y and X modules
        const auto & hit2Geo = trigGeo2.SensitiveVolume(tempHits_F[f_test].channel);
       const auto hit2Center = hit2Geo.GetCenter();
       bool moduleMatched;
       moduleMatched=moduleMatcher(tempHits_F[f_test].module, tempHits_F[f].module);
       if (moduleMatched) {
         // Get the center of the hits (CRT_Res=2.5 cm)
         double hitX = hit1Center.X();
         for (unsigned int a = 0; a < tempHits_F.size(); a++)
         {
         if(tempHits_F[a].module==tempHits_F[f].module && (tempHits_F[a].channel-1)==tempHits_F[f].channel) hitX=hit1Center.X()+1.25;
         }
         double hitYPrelim=hit2Center.Y();
         for (unsigned int a = 0; a < tempHits_F.size(); a++)
         {
         if(tempHits_F[a].module==tempHits_F[f_test].module && (tempHits_F[a].channel-1)==tempHits_F[f_test].channel) hitYPrelim=hit2Center.Y()+1.25;
         }
         
 
         
         double hitY=hitYPrelim;
         double hitZ = (hit1Center.Z() + hit2Center.Z()) / 2.f;
 
         recoHits rHits;
         rHits.hitPositionX = hitX;
         rHits.hitPositionY = hitY;
         rHits.hitPositionZ = hitZ;
         rHits.trigNumberX=tempHits_F[f].triggerNumber;
         rHits.trigNumberY=tempHits_F[f_test].triggerNumber;
         rHits.timeAvg = (tempHits_F[f_test].triggerTime+tempHits_F[f].triggerTime)/2.0;
        primaryHits_F.push_back(rHits); // Add array
     }
     }
   }
   for (unsigned int f = 0; f < tempHits_B.size(); f++) {
     for (unsigned int f_test = 0; f_test < tempHits_B.size(); f_test++) { // Same as above but for back CRT
        if (fabs(tempHits_B[f_test].triggerTime-tempHits_B[f].triggerTime)>fModuletoModuleTimingCut) continue;
 
       const auto & trigGeo = geom -> AuxDet(tempHits_B[f].module);
       const auto & trigGeo2 = geom -> AuxDet(tempHits_B[f_test].module);
       const auto & hit1Geo = trigGeo.SensitiveVolume(tempHits_B[f].channel);
       const auto hit1Center = hit1Geo.GetCenter();
 
       const auto & hit2Geo = trigGeo2.SensitiveVolume(tempHits_B[f_test].channel);
       const auto hit2Center = hit2Geo.GetCenter();
       bool moduleMatched;
       moduleMatched=moduleMatcher(tempHits_B[f_test].module, tempHits_B[f].module);
 
       if (moduleMatched) {
         double hitX = hit1Center.X();
         
         
         for (unsigned int a = 0; a < tempHits_B.size(); a++)
         {
         if(tempHits_B[a].module==tempHits_B[f].module && (tempHits_B[a].channel-1)==tempHits_B[f].channel) hitX=hit1Center.X()+1.25;
         }
         
         double hitYPrelim = hit2Center.Y();
         
         for (unsigned int a = 0; a < tempHits_B.size(); a++)
         {
         if(tempHits_B[a].module==tempHits_B[f_test].module && (tempHits_B[a].channel-1)==tempHits_B[f_test].channel) hitYPrelim=hit2Center.Y()+1.25;
         }
         double hitY=hitYPrelim;
 
         
         double hitZ = (hit1Center.Z() + hit2Center.Z()) / 2.f;
 
         recoHits rHits;
         rHits.hitPositionX = hitX;
         rHits.hitPositionY = hitY;
         rHits.hitPositionZ = hitZ;
         rHits.trigNumberX=tempHits_B[f].triggerNumber;
         rHits.trigNumberY=tempHits_B[f_test].triggerNumber;
         rHits.timeAvg = (tempHits_B[f_test].triggerTime+tempHits_B[f].triggerTime)/2.0;
          primaryHits_B.push_back(rHits); 
 
          }
     }
   }
   vector < art::Ptr < recob::Track > > trackList;
   auto trackListHandle = event.getHandle < vector < recob::Track > >(fTrackModuleLabel);
   if (trackListHandle) {
     art::fill_ptr_vector(trackList, trackListHandle);
   }
   else{
     event.put(std::move(T0col)); event.put(std::move(CRTTrack)); event.put(std::move(TPCCRTassn));  event.put(std::move(TPCT0assn)); event.put(std::move(CRTTriggerassn));
     return;
   }
 
   vector<art::Ptr<recob::PFParticle> > pfplist;
   auto PFPListHandle = event.getHandle< std::vector<recob::PFParticle> >("pandora");
   if (PFPListHandle) art::fill_ptr_vector(pfplist, PFPListHandle);
   if(pfplist.size()<1) return;
   art::FindManyP<anab::T0> trk_t0_assn_v(PFPListHandle, event ,"pandora");
     art::FindManyP<recob::PFParticle> pfp_trk_assn(trackListHandle,event,"pandoraTrack");
   int nTracksReco = trackList.size();
   art::FindManyP < recob::Hit > hitsFromTrack(trackListHandle, event, fTrackModuleLabel);
 
   for (int iRecoTrack = 0; iRecoTrack < nTracksReco; ++iRecoTrack) {
     std::vector< art::Ptr<recob::Hit> > allHits =  hitsFromTrack.at(iRecoTrack);
 
       art::Ptr<recob::Track> ptrack(trackListHandle, iRecoTrack);
      
         std::vector<art::Ptr<recob::PFParticle>> pfps=pfp_trk_assn.at(iRecoTrack);
         if(!pfps.size()) continue;
         std::vector<art::Ptr<anab::T0>> t0s=trk_t0_assn_v.at(pfps[0].key());
     
 
     double trackStartPositionZ_noSCE = trackList[iRecoTrack]->Vertex().Z();
     double trackEndPositionZ_noSCE = trackList[iRecoTrack] -> End().Z();
 
     double trackStartPositionX_notCorrected = trackList[iRecoTrack]->Vertex().X();
     double trackStartPositionY_noSCE = trackList[iRecoTrack]->Vertex().Y();
 
 
     double trackEndPositionX_notCorrected = trackList[iRecoTrack] -> End().X();
     double trackEndPositionY_noSCE = trackList[iRecoTrack] -> End().Y();
 
     int firstHit=0;
     int lastHit=allHits.size()-2;
     if (trackStartPositionZ_noSCE>trackEndPositionZ_noSCE){
     trackEndPositionZ_noSCE = trackList[iRecoTrack]->Vertex().Z();
     trackStartPositionZ_noSCE = trackList[iRecoTrack] -> End().Z();
     trackEndPositionX_notCorrected = trackList[iRecoTrack]->Vertex().X();
     trackEndPositionY_noSCE = trackList[iRecoTrack]->Vertex().Y();
 
 
     trackStartPositionX_notCorrected=trackList[iRecoTrack] -> End().X();
     trackStartPositionY_noSCE = trackList[iRecoTrack] -> End().Y();
     firstHit=lastHit;
     lastHit=0;
     }
 
 
 
    
 
 
 
     if ((trackEndPositionZ_noSCE > 660 && trackStartPositionZ_noSCE < 50) || (trackStartPositionZ_noSCE > 660 && trackEndPositionZ_noSCE < 50)) {
 
       double min_delta = DBL_MAX;
       double best_XF = DBL_MAX;
       double best_YF = DBL_MAX;
       double best_ZF = DBL_MAX;
       double best_XB = DBL_MAX;
       double best_YB = DBL_MAX;
       double best_ZB = DBL_MAX;
       double best_dotProductCos = DBL_MAX;
       double best_deltaXF = DBL_MAX;
       double best_deltaYF = DBL_MAX;
       double best_deltaXB = DBL_MAX;
       double best_deltaYB = DBL_MAX;
       double best_T=DBL_MAX;
       int best_trigXF=0;
       int best_trigYF=0;
       int best_trigXB=0;
       int best_trigYB=0;
 
     for (unsigned int f = 0; f < primaryHits_F.size(); f++) {
         for (unsigned int b = 0; b < primaryHits_B.size(); b++) {
 
       double X1 = primaryHits_F[f].hitPositionX;
       double Y1 = primaryHits_F[f].hitPositionY;
       double Z1 = primaryHits_F[f].hitPositionZ;
       double X2 = primaryHits_B[b].hitPositionX;
       double Y2 = primaryHits_B[b].hitPositionY;
       double Z2= primaryHits_B[b].hitPositionZ;
      
 
                 if (fabs(primaryHits_F[f].timeAvg-primaryHits_B[b].timeAvg)>fFronttoBackTimingCut) continue;
                 double t0=(primaryHits_F[f].timeAvg+primaryHits_B[b].timeAvg)/2.f;
                 int tempId = 0;
                 double xOffset=0;
     
 
                 double trackStartPositionX_noSCE=trackStartPositionX_notCorrected;
                 double trackEndPositionX_noSCE=trackEndPositionX_notCorrected;
                 if (t0s.empty()){
                 int RDOffset=0;
                 if (!fMCCSwitch) RDOffset=111;
                 double ticksOffset=0;
 
                 if (!fMCCSwitch) ticksOffset = (t0+RDOffset)/25.f+detProp.GetXTicksOffset(allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
 
                 else if (fMCCSwitch) ticksOffset = (t0/500.f)+detProp.GetXTicksOffset(allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
                 
                xOffset=detProp.ConvertTicksToX(ticksOffset,allHits[firstHit]->WireID().Plane, allHits[firstHit]->WireID().TPC, allHits[firstHit]->WireID().Cryostat);
                 //double xOffset=.08*ticksOffset
                 
          trackStartPositionX_noSCE=trackStartPositionX_notCorrected-xOffset;
          trackEndPositionX_noSCE=trackEndPositionX_notCorrected-xOffset;
         }
 
    double trackStartPositionX=trackStartPositionX_noSCE;
    double trackStartPositionY=trackStartPositionY_noSCE;
    double trackStartPositionZ=trackStartPositionZ_noSCE;
 
    double trackEndPositionX=trackEndPositionX_noSCE;
    double trackEndPositionY=trackEndPositionY_noSCE;
    double trackEndPositionZ=trackEndPositionZ_noSCE;
    if (fSCECorrection && SCE->EnableCalSpatialSCE()){
 if(geom->PositionToTPCID(geo::Point_t(trackEndPositionX, trackEndPositionY, trackEndPositionZ)).deepestIndex()<13 && geom->PositionToTPCID(geo::Point_t(trackStartPositionX, trackStartPositionY, trackStartPositionZ)).deepestIndex()<13){ 
             auto const & posOffsets_F = SCE->GetCalPosOffsets(geo::Point_t(trackStartPositionX, trackStartPositionY, trackStartPositionZ), geom->PositionToTPCID(geo::Point_t(trackStartPositionX, trackStartPositionY, trackStartPositionZ)).deepestIndex());
             trackStartPositionX -= posOffsets_F.X();
             trackStartPositionY += posOffsets_F.Y();
             trackStartPositionZ += posOffsets_F.Z();
             auto const & posOffsets_B = SCE->GetCalPosOffsets(geo::Point_t(trackEndPositionX, trackEndPositionY, trackEndPositionZ), geom->PositionToTPCID(geo::Point_t(trackEndPositionX, trackEndPositionY, trackEndPositionZ)).deepestIndex());
             trackEndPositionX -= posOffsets_B.X();
             trackEndPositionY += posOffsets_B.Y();
             trackEndPositionZ += posOffsets_B.Z();
         }
     }
 
 
 
 
 
         // Make metrics for a CRT pair to compare later
         TVector3 trackStart(trackStartPositionX, trackStartPositionY, trackStartPositionZ);
         TVector3 trackEnd(trackEndPositionX, trackEndPositionY, trackEndPositionZ);
         TVector3 v1(X1,Y1,Z1);
         TVector3 v2(X2, Y2, Z2);
 
             TVector3 v4(trackStartPositionX,
                         trackStartPositionY,
                         trackStartPositionZ);
             TVector3 v5(trackEndPositionX,
                         trackEndPositionY,
                         trackEndPositionZ);
         TVector3 trackVector = (v5-v4).Unit();
         TVector3 hitVector=(v2-v1).Unit();
 
 
 
 
               double predictedHitPositionY1 = (v1.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.Y()-v5.Y())+v5.Y();
               double predictedHitPositionY2 = (v2.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.Y()-v5.Y())+v5.Y();
 
               double predictedHitPositionX1 = (v1.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.X()-v5.X())+v5.X();
               double predictedHitPositionX2 = (v2.Z()-v5.Z())/(v4.Z()-v5.Z())*(v4.X()-v5.X())+v5.X();
 
         double dotProductCos=trackVector*hitVector;
 
         double deltaX1 = (predictedHitPositionX1-X1);
         double deltaX2 = (predictedHitPositionX2-X2);
 
         double deltaY1 = (predictedHitPositionY1-Y1);
 
         double deltaY2 = (predictedHitPositionY2-Y2);
 
 
 
       tempId++;
 
      //iRecoTrack
       if (min_delta > std::abs(deltaX1)+std::abs(deltaX2) + std::abs(deltaY1)+std::abs(deltaY2) ){
 
            min_delta=std::abs(deltaX1)+std::abs(deltaX2) + std::abs(deltaY1)+std::abs(deltaY2);
             best_XF = X1;
             best_YF = Y1;
             best_ZF = Z1;
             best_XB = X2;
             best_YB = Y2;
             best_ZB = Z2;
             best_dotProductCos = dotProductCos;
             best_deltaXF = deltaX1;
             best_deltaYF = deltaY1;
             best_deltaXB = deltaX2;
             best_deltaYB = deltaY2;
             best_trigXF=primaryHits_F[f].trigNumberX;
             best_trigYF=primaryHits_F[f].trigNumberY;
             best_trigXB=primaryHits_B[b].trigNumberX;
             best_trigYB=primaryHits_B[b].trigNumberY;
             best_T = t0;
             if (!fMCCSwitch) best_T=(111.f+best_T)*20.f;
             // Added 111 tick CRT-CTB offset
           }
         }
       }
       if (std::abs(best_dotProductCos)>0.99 && std::abs(best_deltaXF)+std::abs(best_deltaXB)<40 && std::abs(best_deltaYF)+std::abs(best_deltaYB)<40 ) {
         //std::cout<<"Found match with TPC*CRT "<<best_dotProductCos<<std::endl;
 
 //std::cout<<"Displacement in front and back "<<best_deltaXF<<","<<best_deltaYF<<","<<best_deltaXB<<","<<best_deltaYB<<std::endl;
         std::vector<float> hitF;
         std::vector<float> hitB;
         hitF.push_back(best_XF); hitF.push_back(best_YF); hitF.push_back(best_ZF);
         hitB.push_back(best_XB); hitB.push_back(best_YB); hitB.push_back(best_ZB);
 
         CRTTrack->push_back(anab::CosmicTag(hitF,hitB, std::abs(best_dotProductCos),anab::CosmicTagID_t::kNotTagged));
 
         T0col->push_back(anab::T0(best_T, 13,2,iRecoTrack,best_dotProductCos));
         util::CreateAssn(*this, event, *CRTTrack, trackList[iRecoTrack], *TPCCRTassn);
         util::CreateAssn(*this, event, *T0col, trackList[iRecoTrack], *TPCT0assn);
         util::CreateAssn(*this, event, *CRTTrack, crtList[best_trigXF], *CRTTriggerassn);
         util::CreateAssn(*this, event, *CRTTrack, crtList[best_trigYF], *CRTTriggerassn);
 
         util::CreateAssn(*this, event, *CRTTrack, crtList[best_trigXB], *CRTTriggerassn);
 
         util::CreateAssn(*this, event, *CRTTrack, crtList[best_trigYB], *CRTTriggerassn);
         
         }
       }
     }
 
  event.put(std::move(T0col)); event.put(std::move(CRTTrack)); event.put(std::move(TPCCRTassn));  event.put(std::move(TPCT0assn)); event.put(std::move(CRTTriggerassn));
 }