LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CosmicTrackTagger_module.cc
Go to the documentation of this file.
1 // Class: CosmicTagger
3 // Module Type: producer
4 // File: CosmicTrackTagger_module.cc
5 //
6 // Generated at Mon Sep 24 18:21:00 2012 by Sarah Lockwitz using artmod
7 // from art v1_02_02.
8 // artmod -e beginJob -e reconfigure -e endJob producer trkf::CosmicTrackTagger
10 
16 
17 #include <iostream>
18 
23 
28 
29 #include "TStopwatch.h"
30 
31 namespace cosmic {
32  class CosmicTrackTagger;
33 }
34 
36 public:
37  explicit CosmicTrackTagger(fhicl::ParameterSet const& p);
38 
39  void produce(art::Event& e) override;
40 
41 private:
42  std::string fTrackModuleLabel;
48 };
49 
51 {
52  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
53  auto const detp =
55  auto const* geo = lar::providerFrom<geo::Geometry>();
56 
57  fDetHalfHeight = geo->DetHalfHeight();
58  fDetWidth = 2. * geo->DetHalfWidth();
59  fDetLength = geo->DetLength();
60 
61  float fSamplingRate = sampling_rate(clock_data);
62 
63  fTrackModuleLabel = p.get<std::string>("TrackModuleLabel", "track");
64  fEndTickPadding = p.get<int>("EndTickPadding", 50);
65 
66  fTPCXBoundary = p.get<float>("TPCXBoundary", 5);
67  fTPCYBoundary = p.get<float>("TPCYBoundary", 5);
68  fTPCZBoundary = p.get<float>("TPCZBoundary", 5);
69 
70  const double driftVelocity = detp.DriftVelocity(detp.Efield(), detp.Temperature()); // cm/us
71 
73  2 * geo->DetHalfWidth() / (driftVelocity * fSamplingRate / 1000); // ~3200 for uB
74  fMinTickDrift = clock_data.Time2Tick(clock_data.TriggerTime());
76 
77  produces<std::vector<anab::CosmicTag>>();
78  produces<art::Assns<recob::Track, anab::CosmicTag>>();
79 }
80 
82 {
83  // Implementation of required member function here.
84 
85  std::unique_ptr<std::vector<anab::CosmicTag>> cosmicTagTrackVector(
86  new std::vector<anab::CosmicTag>);
87  std::unique_ptr<art::Assns<recob::Track, anab::CosmicTag>> assnOutCosmicTagTrack(
89 
90  TStopwatch ts;
91 
93  e.getByLabel(fTrackModuleLabel, Trk_h);
94  std::vector<art::Ptr<recob::Track>> TrkVec;
95  art::fill_ptr_vector(TrkVec, Trk_h);
96 
98  // LOOPING OVER INSPILL TRACKS
100 
101  art::FindManyP<recob::Hit> hitsSpill(Trk_h, e, fTrackModuleLabel);
102 
103  for (unsigned int iTrack = 0; iTrack < Trk_h->size(); iTrack++) {
104 
105  int isCosmic = 0;
107 
108  art::Ptr<recob::Track> tTrack = TrkVec.at(iTrack);
109  std::vector<art::Ptr<recob::Hit>> HitVec = hitsSpill.at(iTrack);
110 
111  if (iTrack != tTrack.key()) { std::cout << "Mismatch in track index/key" << std::endl; }
112 
113  // A BETTER WAY OF FINDING END POINTS:
114  auto tVector1 = tTrack->Vertex();
115  auto tVector2 = tTrack->End();
116 
117  float trackEndPt1_X = tVector1.X();
118  float trackEndPt1_Y = tVector1.Y();
119  float trackEndPt1_Z = tVector1.Z();
120  float trackEndPt2_X = tVector2.X();
121  float trackEndPt2_Y = tVector2.Y();
122  float trackEndPt2_Z = tVector2.Z();
123 
124  if (trackEndPt1_X != trackEndPt1_X || trackEndPt1_Y != trackEndPt1_Y ||
125  trackEndPt1_Z != trackEndPt1_Z || trackEndPt2_X != trackEndPt2_X ||
126  trackEndPt2_Y != trackEndPt2_Y || trackEndPt2_Z != trackEndPt2_Z) {
127  std::cerr << "!!! FOUND A PROBLEM... the length is: " << tTrack->Length()
128  << " np: " << tTrack->NumberTrajectoryPoints() << " id: " << tTrack->ID() << " "
129  << tTrack << std::endl;
130  std::vector<float> tempPt1, tempPt2;
131  tempPt1.push_back(-999);
132  tempPt1.push_back(-999);
133  tempPt1.push_back(-999);
134  tempPt2.push_back(-999);
135  tempPt2.push_back(-999);
136  tempPt2.push_back(-999);
137  cosmicTagTrackVector->emplace_back(tempPt1, tempPt2, -999, tag_id);
138  util::CreateAssn(e, *cosmicTagTrackVector, tTrack, *assnOutCosmicTagTrack);
139  continue; // I don't want to deal with these "tracks"
140  }
141 
143  // Getting first and last ticks
145  float tick1 = 9999;
146  float tick2 = -9999;
147  bool dumpMe(false);
148 
149  for (unsigned int p = 0; p < HitVec.size(); p++) {
150  if (dumpMe) {
151  std::cout << "###>> Hit key: " << HitVec[p].key()
152  << ", peak - RMS: " << HitVec[p]->PeakTimeMinusRMS()
153  << ", peak + RMS: " << HitVec[p]->PeakTimePlusRMS() << std::endl;
154  }
155  if (HitVec[p]->PeakTimeMinusRMS() < tick1) tick1 = HitVec[p]->PeakTimeMinusRMS();
156  if (HitVec[p]->PeakTimePlusRMS() > tick2) tick2 = HitVec[p]->PeakTimePlusRMS();
157  }
158 
160  // Are any of the ticks outside of the ReadOutWindow ?
162  if (tick1 < fMinTickDrift || tick2 > fMaxTickDrift) {
163  isCosmic = 1;
165  }
166 
168  // Now check Y & Z boundaries:
170  int nBdY = 0, nBdZ = 0;
171  if (isCosmic == 0) {
172 
173  // Checking lower side of TPC
174  if (fabs(fDetHalfHeight + trackEndPt1_Y) < fTPCYBoundary ||
175  fabs(fDetHalfHeight + trackEndPt2_Y) < fTPCYBoundary || trackEndPt1_Y < -fDetHalfHeight ||
176  trackEndPt2_Y < -fDetHalfHeight)
177  nBdY++;
178 
179  // Checking upper side of TPC
180  if (fabs(fDetHalfHeight - trackEndPt1_Y) < fTPCYBoundary ||
181  fabs(fDetHalfHeight - trackEndPt2_Y) < fTPCYBoundary || trackEndPt1_Y > fDetHalfHeight ||
182  trackEndPt2_Y > fDetHalfHeight)
183  nBdY++;
184 
185  if (fabs(trackEndPt1_Z - fDetLength) < fTPCZBoundary ||
186  fabs(trackEndPt2_Z - fDetLength) < fTPCZBoundary)
187  nBdZ++;
188  if (fabs(trackEndPt1_Z) < fTPCZBoundary || fabs(trackEndPt2_Z) < fTPCZBoundary) nBdZ++;
189  if ((nBdY + nBdZ) > 1) {
190  isCosmic = 2;
191  if (nBdY > 1)
193  else if (nBdZ > 1)
195  else
197  }
198  else if ((nBdY + nBdZ) == 1) {
199  isCosmic = 3;
200  if (nBdY == 1)
202  else if (nBdZ == 1)
204  }
205  }
206 
207  std::vector<float> endPt1;
208  std::vector<float> endPt2;
209  endPt1.push_back(trackEndPt1_X);
210  endPt1.push_back(trackEndPt1_Y);
211  endPt1.push_back(trackEndPt1_Z);
212  endPt2.push_back(trackEndPt2_X);
213  endPt2.push_back(trackEndPt2_Y);
214  endPt2.push_back(trackEndPt2_Z);
215 
216  float cosmicScore = isCosmic > 0 ? 1 : 0;
217  if (isCosmic == 3) cosmicScore = 0.5;
218 
220  // Doing a very basic check on X boundaries
221  // this gets the types of tracks that go through both X boundaries of the detector
222  if (fabs(trackEndPt1_X - trackEndPt2_X) > fDetWidth - fTPCXBoundary) {
223  cosmicScore = 1;
224  isCosmic = 4;
226  }
227 
228  cosmicTagTrackVector->emplace_back(endPt1, endPt2, cosmicScore, tag_id);
229 
230  util::CreateAssn(e, *cosmicTagTrackVector, tTrack, *assnOutCosmicTagTrack);
231  }
232  // END OF LOOPING OVER INSPILL TRACKS
233 
237  float dE = 0, dS = 0, temp = 0, IScore = 0;
238  unsigned int IndexE = 0, iTrk1 = 0, iTrk = 0;
240 
241  for (iTrk = 0; iTrk < Trk_h->size(); iTrk++) {
242  art::Ptr<recob::Track> tTrk = TrkVec.at(iTrk);
243  if ((*cosmicTagTrackVector)[iTrk].CosmicScore() == 0) {
244  auto tStart = tTrk->Vertex();
245  auto tEnd = tTrk->End();
246  unsigned int l = 0;
247  for (iTrk1 = 0; iTrk1 < Trk_h->size(); iTrk1++) {
248  art::Ptr<recob::Track> tTrk1 = TrkVec.at(iTrk1);
249  float getScore = (*cosmicTagTrackVector)[iTrk1].CosmicScore();
250  if (getScore == 1 || getScore == 0.5) {
251  anab::CosmicTagID_t getType = (*cosmicTagTrackVector)[iTrk1].CosmicType();
252  auto tStart1 = tTrk1->Vertex();
253  auto tEnd1 = tTrk1->End();
254  auto NumE = (tEnd - tStart1).Cross(tEnd - tEnd1);
255  auto DenE = tEnd1 - tStart1;
256  dE = NumE.R() / DenE.R();
257  if (l == 0) {
258  temp = dE;
259  IndexE = iTrk1;
260  IScore = getScore;
261  IType = getType;
262  }
263  if (dE < temp) {
264  temp = dE;
265  IndexE = iTrk1;
266  IScore = getScore;
267  IType = getType;
268  }
269  l++;
270  }
271  } //End Trk1 loop
272  art::Ptr<recob::Track> tTrkI = TrkVec.at(IndexE);
273  auto tStartI = tTrkI->Vertex();
274  auto tEndI = tTrkI->End();
275  auto NumS = (tStart - tStartI).Cross(tStart - tEndI);
276  auto DenS = tEndI - tStartI;
277  dS = NumS.R() / DenS.R();
278  if (((dS < 5 && temp < 5) || (dS < temp && dS < 5)) && (tTrk->Length() < 60)) {
279  (*cosmicTagTrackVector)[iTrk].CosmicScore() = IScore - 0.05;
280  (*cosmicTagTrackVector)[iTrk].CosmicType() = IType;
281  }
282  } // end cosmicScore==0 loop
283  } // end iTrk loop
284 
285  e.put(std::move(cosmicTagTrackVector));
286  e.put(std::move(assnOutCosmicTagTrack));
287 
288  TrkVec.clear();
289 
290 } // end of produce
291 
std::string getType(cet::LibraryManager const &lm, std::string const &fullSpec)
Utilities related to art service access.
enum anab::cosmic_tag_id CosmicTagID_t
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:136
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:207
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
key_type key() const noexcept
Definition: Ptr.h:166
Point_t const & Vertex() const
Access to track position at different points.
Definition: Track.h:158
Provides recob::Track data product.
void produce(art::Event &e) override
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
int ID() const
Definition: Track.h:244
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
Point_t const & End() const
Access to track position at different points.
Definition: Track.h:159
CosmicTrackTagger(fhicl::ParameterSet const &p)
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
Float_t e
Definition: plot.C:35
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Namespace collecting geometry-related classes utilities.
art framework interface to geometry description