All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes | Public Types | Public Member Functions | Private Attributes | List of all members
sbnd::CRTT0MatchingAna Class Reference
Inheritance diagram for sbnd::CRTT0MatchingAna:

Classes

struct  Config
 

Public Types

using Parameters = art::EDAnalyzer::Table< Config >
 

Public Member Functions

 CRTT0MatchingAna (Parameters const &config)
 
virtual void beginJob () override
 
virtual void analyze (const art::Event &event) override
 
virtual void endJob () override
 
double DistToCrtHit (TVector3 trackPos, sbn::crt::CRTHit crtHit)
 

Private Attributes

art::InputTag fSimModuleLabel
 name of detsim producer More...
 
art::InputTag fCRTHitLabel
 name of CRT producer More...
 
art::InputTag fTPCTrackLabel
 name of CRT producer More...
 
bool fVerbose
 print information about what's going on More...
 
CRTT0MatchAlg t0Alg
 
CRTGeoAlg fCrtGeo
 
TPCGeoAlg fTpcGeo
 
CRTBackTracker fCrtBackTrack
 
std::map< std::string, TH1D * > hDCA
 
std::map< std::string, TH1D * > hMatchDCA
 
std::map< std::string, TH1D * > hNoMatchDCA
 
std::map< std::string, TH1D * > hDoL
 
std::map< std::string, TH1D * > hMatchDoL
 
std::map< std::string, TH1D * > hNoMatchDoL
 
std::map< std::string, TH1D * > hT0
 
std::map< std::string, TH1D * > hMatchT0
 
std::map< std::string, TH1D * > hNoMatchT0
 
std::map< std::string, TH1D * > hEffDCATotal
 
std::map< std::string, TH1D * > hEffDCAReco
 
std::map< std::string, TH1D * > hEffDoLTotal
 
std::map< std::string, TH1D * > hEffDoLReco
 
std::map< std::string, TH1D * > hEffLengthTotal
 
std::map< std::string, TH1D * > hEffLengthReco
 
std::map< std::string, TH1D * > hPurityDCATotal
 
std::map< std::string, TH1D * > hPurityDCAReco
 
std::map< std::string, TH1D * > hPurityDoLTotal
 
std::map< std::string, TH1D * > hPurityDoLReco
 
std::map< std::string, TH1D * > hPurityLengthTotal
 
std::map< std::string, TH1D * > hPurityLengthReco
 

Detailed Description

Definition at line 62 of file sbndcode/sbndcode/CRT/CRTAna/CRTT0MatchingAna_module.cc.

Member Typedef Documentation

using sbnd::CRTT0MatchingAna::Parameters = art::EDAnalyzer::Table<Config>

Constructor & Destructor Documentation

sbnd::CRTT0MatchingAna::CRTT0MatchingAna ( Parameters const &  config)
explicit

Definition at line 164 of file sbndcode/sbndcode/CRT/CRTAna/CRTT0MatchingAna_module.cc.

165  : EDAnalyzer(config)
166  , fSimModuleLabel (config().SimModuleLabel())
167  , fCRTHitLabel (config().CRTHitLabel())
168  , fTPCTrackLabel (config().TPCTrackLabel())
169  , fVerbose (config().Verbose())
170  , t0Alg (config().CRTT0Alg())
171  , fCrtBackTrack (config().CrtBackTrack())
172  {
173 
174  } //CRTT0MatchingAna()
BEGIN_PROLOG Verbose
BEGIN_PROLOG don t mess with this TPCTrackLabel
bool fVerbose
print information about what&#39;s going on

Member Function Documentation

void sbnd::CRTT0MatchingAna::analyze ( const art::Event &  event)
overridevirtual

Definition at line 220 of file sbndcode/sbndcode/CRT/CRTAna/CRTT0MatchingAna_module.cc.

221  {
222 
223  // Fetch basic event info
224  if(fVerbose){
225  std::cout<<"============================================"<<std::endl
226  <<"Run = "<<event.run()<<", SubRun = "<<event.subRun()<<", Event = "<<event.id().event()<<std::endl
227  <<"============================================"<<std::endl;
228  }
229 
230  //----------------------------------------------------------------------------------------------------------
231  // GETTING PRODUCTS
232  //----------------------------------------------------------------------------------------------------------
233 
234  // Get g4 particles
235  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
236  auto particleHandle = event.getValidHandle<std::vector<simb::MCParticle>>(fSimModuleLabel);
237 
238  // Get CRT hits from the event
239  art::Handle< std::vector<sbn::crt::CRTHit>> crtHitHandle;
240  std::vector<art::Ptr<sbn::crt::CRTHit> > crtHitList;
241  if (event.getByLabel(fCRTHitLabel, crtHitHandle))
242  art::fill_ptr_vector(crtHitList, crtHitHandle);
243 
244  // Get reconstructed tracks from the event
245  auto tpcTrackHandle = event.getValidHandle<std::vector<recob::Track>>(fTPCTrackLabel);
246  art::FindManyP<recob::Hit> findManyHits(tpcTrackHandle, event, fTPCTrackLabel);
247 
248  fCrtBackTrack.Initialize(event);
249 
250  //----------------------------------------------------------------------------------------------------------
251  // TRUTH MATCHING
252  //----------------------------------------------------------------------------------------------------------
253 
254  std::map<int, simb::MCParticle> particles;
255  // Loop over the true particles
256  for (auto const& particle: (*particleHandle)){
257 
258  // Make map with ID
259  int partID = particle.TrackId();
260  particles[partID] = particle;
261 
262  }
263 
264  std::map<int, std::vector<std::string>> crtTaggerMap;
265  std::vector<sbn::crt::CRTHit> crtHits;
266  int hit_i = 0;
267  double minHitTime = 99999;
268  double maxHitTime = -99999;
269 
270  for(auto const& hit : (*crtHitHandle)){
271  double hitTime = (double)(int)hit.ts1_ns * 1e-3;
272  if(hitTime < minHitTime) minHitTime = hitTime;
273  if(hitTime > maxHitTime) maxHitTime = hitTime;
274 
275  crtHits.push_back(hit);
276  int trueId = fCrtBackTrack.TrueIdFromHitId(event, hit_i);
277  hit_i++;
278  if(trueId == -99999) continue;
279  if(crtTaggerMap.find(trueId) == crtTaggerMap.end()){
280  crtTaggerMap[trueId].push_back(hit.tagger);
281  }
282  else if(std::find(crtTaggerMap[trueId].begin(), crtTaggerMap[trueId].end(), hit.tagger) == crtTaggerMap[trueId].end()){
283  crtTaggerMap[trueId].push_back(hit.tagger);
284  }
285  }
286 
287  // std::cout << " New event" << std::endl;
288  //----------------------------------------------------------------------------------------------------------
289  // DISTANCE OF CLOSEST APPROACH ANALYSIS
290  //----------------------------------------------------------------------------------------------------------
291 
292  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
293  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(event, clockData);
294 
295  // Loop over reconstructed tracks
296  for (auto const& tpcTrack : (*tpcTrackHandle)){
297  // Get the associated hits
298  std::vector<art::Ptr<recob::Hit>> hits = findManyHits.at(tpcTrack.ID());
299  int trackTrueID = RecoUtils::TrueParticleIDFromTotalRecoHits(clockData, hits, false);
300  if(particles.find(trackTrueID) == particles.end()) continue;
301  // Only consider primary muons
302  if(!(std::abs(particles[trackTrueID].PdgCode()) == 13 && particles[trackTrueID].Mother() == 0)) continue;
303 
304  // Only consider particles inside the reco hit time window
305  double trueTime = particles[trackTrueID].T() * 1e-3;
306  if(trueTime < minHitTime || trueTime > maxHitTime) continue;
307 
308  // std::cout << "new track " << trueTime << std::endl;
309  // Calculate t0 from CRT Hit matching
310  matchCand closest = t0Alg.GetClosestCRTHit(detProp, tpcTrack, crtHits, event);
311  // std::cout << "closest match " << closest.t0 << std::endl;
312  double sin_angle = -99999;
313  if(closest.dca != -99999){
314  hDCA[closest.thishit.tagger]->Fill(closest.dca);
315  hDCA["All"]->Fill(closest.dca);
316  sin_angle = closest.dca/closest.extrapLen;
317  hDoL[closest.thishit.tagger]->Fill(sin_angle);
318  hDoL["All"]->Fill(sin_angle);
319  hT0[closest.thishit.tagger]->Fill(closest.t0);
320  hT0["All"]->Fill(closest.t0);
321 
322 
323  // Is hit matched to that track
324  int hitTrueID = fCrtBackTrack.TrueIdFromTotalEnergy(event, closest.thishit);
325  if(hitTrueID == trackTrueID && hitTrueID != -99999){
326  hMatchDCA[closest.thishit.tagger]->Fill(closest.dca);
327  hMatchDCA["All"]->Fill(closest.dca);
328  hMatchDoL[closest.thishit.tagger]->Fill(sin_angle);
329  hMatchDoL["All"]->Fill(sin_angle);
330  hMatchT0[closest.thishit.tagger]->Fill(closest.t0);
331  hMatchT0["All"]->Fill(closest.t0);
332  }
333  else{
334  hNoMatchDCA[closest.thishit.tagger]->Fill(closest.dca);
335  hNoMatchDCA["All"]->Fill(closest.dca);
336  hNoMatchDoL[closest.thishit.tagger]->Fill(sin_angle);
337  hNoMatchDoL["All"]->Fill(sin_angle);
338  hNoMatchT0[closest.thishit.tagger]->Fill(closest.t0);
339  hNoMatchT0["All"]->Fill(closest.t0);
340  }
341  }
342 
343  int nbins = hEffDCATotal.begin()->second->GetNbinsX();
344  for(int i = 0; i < nbins; i++){
345  double DCAcut = hEffDCATotal.begin()->second->GetBinCenter(i);
346 
347  // Fill total efficiency histogram with each cut if track matches any hits
348  if(crtTaggerMap.find(trackTrueID) != crtTaggerMap.end()){
349  for(auto const& tagger : crtTaggerMap[trackTrueID]){
350  hEffDCATotal[tagger]->Fill(DCAcut);
351 
352  // If closest hit is below limit and track matches any hits then fill efficiency
353  if(closest.dca < DCAcut && closest.dca != -99999){
354  hEffDCAReco[tagger]->Fill(DCAcut);
355  }
356  }
357  // Fill total efficiency histograms
358  hEffDCATotal["All"]->Fill(DCAcut);
359  if(closest.dca < DCAcut && closest.dca != -99999){
360  hEffDCAReco["All"]->Fill(DCAcut);
361  }
362  }
363 
364  // Fill total purity histogram with each cut if closest hit is below limit
365  if(closest.dca < DCAcut && closest.dca != -99999){
366  hPurityDCATotal[closest.thishit.tagger]->Fill(DCAcut);
367  hPurityDCATotal["All"]->Fill(DCAcut);
368 
369  // If closest hit is below limit and matched time is correct then fill purity
370  double hitTime = closest.thishit.ts1_ns * 1e-3;
371  if(particles.find(trackTrueID) != particles.end()){
372  if(std::abs(hitTime - trueTime) < 2.){
373  hPurityDCAReco[closest.thishit.tagger]->Fill(DCAcut);
374  hPurityDCAReco["All"]->Fill(DCAcut);
375  }
376  }
377 
378  }
379  }
380 
381  nbins = hEffDoLTotal.begin()->second->GetNbinsX();
382 
383  for(int i = 0; i < nbins; i++){
384  double DCAcut = hEffDoLTotal.begin()->second->GetBinCenter(i);
385 
386  // Fill total efficiency histogram with each cut if track matches any hits
387  if(crtTaggerMap.find(trackTrueID) != crtTaggerMap.end()){
388  for(auto const& tagger : crtTaggerMap[trackTrueID]){
389  hEffDoLTotal[tagger]->Fill(DCAcut);
390 
391  // If closest hit is below limit and track matches any hits then fill efficiency
392  if(sin_angle < DCAcut && closest.dca != -99999){
393  hEffDoLReco[tagger]->Fill(DCAcut);
394  }
395  }
396  // Fill total efficiency histograms
397  hEffDoLTotal["All"]->Fill(DCAcut);
398  if(sin_angle < DCAcut && closest.dca != -99999){
399  hEffDoLReco["All"]->Fill(DCAcut);
400  }
401  }
402 
403  // Fill total purity histogram with each cut if closest hit is below limit
404  if(sin_angle < DCAcut && closest.dca != -99999){
405  hPurityDoLTotal[closest.thishit.tagger]->Fill(DCAcut);
406  hPurityDoLTotal["All"]->Fill(DCAcut);
407 
408  // If closest hit is below limit and matched time is correct then fill purity
409  double hitTime = closest.thishit.ts1_ns * 1e-3;
410  if(particles.find(trackTrueID) != particles.end()){
411  if(std::abs(hitTime - trueTime) < 2.){
412  hPurityDoLReco[closest.thishit.tagger]->Fill(DCAcut);
413  hPurityDoLReco["All"]->Fill(DCAcut);
414  }
415  }
416 
417  }
418  }
419 
420  double fixedCut = 30.;
421 
422  // Fill total efficiency histogram with each cut if track matches any hits
423  if(crtTaggerMap.find(trackTrueID) != crtTaggerMap.end()){
424  for(auto const& tagger : crtTaggerMap[trackTrueID]){
425  hEffLengthTotal[tagger]->Fill(tpcTrack.Length());
426 
427  // If closest hit is below limit and track matches any hits then fill efficiency
428  if(closest.dca < fixedCut && closest.dca >=0 ){
429  hEffLengthReco[tagger]->Fill(tpcTrack.Length());
430  }
431  }
432  // Fill total efficiency histograms
433  hEffLengthTotal["All"]->Fill(tpcTrack.Length());
434  if(closest.dca < fixedCut && closest.dca >=0){
435  hEffLengthReco["All"]->Fill(tpcTrack.Length());
436  }
437  }
438 
439  // Fill total purity histogram with each cut if closest hit is below limit
440  if(closest.dca < fixedCut && closest.dca >= 0){
441  hPurityLengthTotal[closest.thishit.tagger]->Fill(tpcTrack.Length());
442  hPurityLengthTotal["All"]->Fill(tpcTrack.Length());
443 
444  // If closest hit is below limit and matched time is correct then fill purity
445  double hitTime = closest.thishit.ts1_ns * 1e-3;
446  if(particles.find(trackTrueID) != particles.end()){
447  double trueTime = particles[trackTrueID].T() * 1e-3;
448  if(std::abs(hitTime - trueTime) < 2.){
449  hPurityLengthReco[closest.thishit.tagger]->Fill(tpcTrack.Length());
450  hPurityLengthReco["All"]->Fill(tpcTrack.Length());
451  }
452  }
453 
454  }
455 
456  }
457 
458 
459  } // CRTT0MatchingAna::analyze()
matchCand GetClosestCRTHit(detinfo::DetectorPropertiesData const &detProp, recob::Track tpcTrack, std::pair< double, double > t0MinMax, std::vector< sbn::crt::CRTHit > crtHits, int driftDirection)
process_name hit
Definition: cheaterreco.fcl:51
int TrueIdFromTotalEnergy(const art::Event &event, const sbnd::crt::CRTData &data)
int TrueIdFromHitId(const art::Event &event, int hit_i)
T abs(T value)
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
int TrueParticleIDFromTotalRecoHits(detinfo::DetectorClocksData const &clockData, const std::vector< art::Ptr< recob::Hit > > &hits, bool rollup_unsaved_ids=1)
auto begin(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:573
bool fVerbose
print information about what&#39;s going on
do i e
BEGIN_PROLOG could also be cout
auto const detProp
void sbnd::CRTT0MatchingAna::beginJob ( )
overridevirtual

Definition at line 177 of file sbndcode/sbndcode/CRT/CRTAna/CRTT0MatchingAna_module.cc.

178  {
179 
180  // Access tfileservice to handle creating and writing histograms
181  art::ServiceHandle<art::TFileService> tfs;
182  for(size_t i = 0; i < fCrtGeo.NumTaggers() + 1; i++){
183  std::string tagger = "All";
184  if(i < fCrtGeo.NumTaggers()){
185  tagger = fCrtGeo.GetTagger(i).name;
186  }
187  hDCA[tagger] = tfs->make<TH1D>(Form("DCA_%s", tagger.c_str()), "", 50, 0, 100);
188  hMatchDCA[tagger] = tfs->make<TH1D>(Form("MatchDCA_%s", tagger.c_str()), "", 50, 0, 100);
189  hNoMatchDCA[tagger] = tfs->make<TH1D>(Form("NoMatchDCA_%s", tagger.c_str()), "", 50, 0, 100);
190 
191  hDoL[tagger] = tfs->make<TH1D>(Form("DoL_%s", tagger.c_str()), "", 100, 0, 0.25);
192  hMatchDoL[tagger] = tfs->make<TH1D>(Form("MatchDoL_%s", tagger.c_str()), "", 100, 0, 0.25);
193  hNoMatchDoL[tagger] = tfs->make<TH1D>(Form("NoMatchDoL_%s", tagger.c_str()), "", 100, 0, 0.25);
194 
195  hT0[tagger] = tfs->make<TH1D>(Form("T0_%s", tagger.c_str()), "", 600, -3000, 3000);
196  hMatchT0[tagger] = tfs->make<TH1D>(Form("MatchT0_%s", tagger.c_str()), "", 600, -3000, 3000);
197  hNoMatchT0[tagger] = tfs->make<TH1D>(Form("NoMatchT0_%s", tagger.c_str()), "", 600, -3000, 3000);
198 
199  hEffDCATotal[tagger] = tfs->make<TH1D>(Form("EffDCATotal_%s", tagger.c_str()), "", 50, 0, 100);
200  hEffDCAReco[tagger] = tfs->make<TH1D>(Form("EffDCAReco_%s", tagger.c_str()), "", 50, 0, 100);
201  hEffDoLTotal[tagger] = tfs->make<TH1D>(Form("EffDoLTotal_%s", tagger.c_str()), "", 100, 0, 0.25);
202  hEffDoLReco[tagger] = tfs->make<TH1D>(Form("EffDoLReco_%s", tagger.c_str()), "", 100, 0, 0.25);
203  hEffLengthTotal[tagger] = tfs->make<TH1D>(Form("EffLengthTotal_%s", tagger.c_str()), "", 20, 0, 600);
204  hEffLengthReco[tagger] = tfs->make<TH1D>(Form("EffLengthReco_%s", tagger.c_str()), "", 20, 0, 600);
205 
206  hPurityDCATotal[tagger] = tfs->make<TH1D>(Form("PurityDCATotal_%s", tagger.c_str()), "", 50, 0, 100);
207  hPurityDCAReco[tagger] = tfs->make<TH1D>(Form("PurityDCAReco_%s", tagger.c_str()), "", 50, 0, 100);
208  hPurityDoLTotal[tagger] = tfs->make<TH1D>(Form("PurityDoLTotal_%s", tagger.c_str()), "", 100, 0, 0.25);
209  hPurityDoLReco[tagger] = tfs->make<TH1D>(Form("PurityDoLReco_%s", tagger.c_str()), "", 100, 0, 0.25);
210  hPurityLengthTotal[tagger] = tfs->make<TH1D>(Form("PurityLengthTotal_%s", tagger.c_str()), "", 20, 0, 600);
211  hPurityLengthReco[tagger] = tfs->make<TH1D>(Form("PurityLengthReco_%s", tagger.c_str()), "", 20, 0, 600);
212  }
213 
214  // Initial output
215  if(fVerbose) std::cout<<"----------------- CRT T0 Matching Ana Module -------------------"<<std::endl;
216 
217  } // CRTT0MatchingAna::beginJob()
CRTTaggerGeo GetTagger(std::string taggerName) const
bool fVerbose
print information about what&#39;s going on
art::ServiceHandle< art::TFileService > tfs
BEGIN_PROLOG could also be cout
double sbnd::CRTT0MatchingAna::DistToCrtHit ( TVector3  trackPos,
sbn::crt::CRTHit  crtHit 
)
void sbnd::CRTT0MatchingAna::endJob ( )
overridevirtual

Definition at line 462 of file sbndcode/sbndcode/CRT/CRTAna/CRTT0MatchingAna_module.cc.

462  {
463 
464 
465  } // CRTT0MatchingAna::endJob()

Member Data Documentation

CRTBackTracker sbnd::CRTT0MatchingAna::fCrtBackTrack
private
CRTGeoAlg sbnd::CRTT0MatchingAna::fCrtGeo
private
art::InputTag sbnd::CRTT0MatchingAna::fCRTHitLabel
private

name of CRT producer

Definition at line 122 of file sbndcode/sbndcode/CRT/CRTAna/CRTT0MatchingAna_module.cc.

art::InputTag sbnd::CRTT0MatchingAna::fSimModuleLabel
private

name of detsim producer

Definition at line 121 of file sbndcode/sbndcode/CRT/CRTAna/CRTT0MatchingAna_module.cc.

TPCGeoAlg sbnd::CRTT0MatchingAna::fTpcGeo
private
art::InputTag sbnd::CRTT0MatchingAna::fTPCTrackLabel
private

name of CRT producer

Definition at line 123 of file sbndcode/sbndcode/CRT/CRTAna/CRTT0MatchingAna_module.cc.

bool sbnd::CRTT0MatchingAna::fVerbose
private

print information about what's going on

Definition at line 124 of file sbndcode/sbndcode/CRT/CRTAna/CRTT0MatchingAna_module.cc.

std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hDCA
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hDoL
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hEffDCAReco
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hEffDCATotal
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hEffDoLReco
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hEffDoLTotal
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hEffLengthReco
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hEffLengthTotal
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hMatchDCA
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hMatchDoL
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hMatchT0
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hNoMatchDCA
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hNoMatchDoL
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hNoMatchT0
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hPurityDCAReco
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hPurityDCATotal
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hPurityDoLReco
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hPurityDoLTotal
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hPurityLengthReco
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hPurityLengthTotal
private
std::map<std::string, TH1D*> sbnd::CRTT0MatchingAna::hT0
private
CRTT0MatchAlg sbnd::CRTT0MatchingAna::t0Alg
private

The documentation for this class was generated from the following file: