Public Member Functions | Private Member Functions | Private Attributes | List of all members
T0RecoSCECalibrations Class Reference
Inheritance diagram for T0RecoSCECalibrations:
art::EDAnalyzer art::detail::Analyzer art::detail::LegacyModule art::Observer art::ModuleBase

Public Member Functions

 T0RecoSCECalibrations (fhicl::ParameterSet const &p)
 
 T0RecoSCECalibrations (T0RecoSCECalibrations const &)=delete
 
 T0RecoSCECalibrations (T0RecoSCECalibrations &&)=delete
 
T0RecoSCECalibrationsoperator= (T0RecoSCECalibrations const &)=delete
 
T0RecoSCECalibrationsoperator= (T0RecoSCECalibrations &&)=delete
 
void beginJob () override
 
void analyze (art::Event const &e) override
 
- Public Member Functions inherited from art::EDAnalyzer
 EDAnalyzer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDAnalyzer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Analyzer
virtual ~Analyzer () noexcept
 
 Analyzer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 Analyzer (Table< Config > const &config)
 
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::Observer
 ~Observer () noexcept
 
 Observer (Observer const &)=delete
 
 Observer (Observer &&)=delete
 
Observeroperator= (Observer const &)=delete
 
Observeroperator= (Observer &&)=delete
 
void registerProducts (ProductDescriptions &, ModuleDescription const &)
 
void fillDescriptions (ModuleDescription const &)
 
fhicl::ParameterSetID selectorConfig () const
 
- 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)
 

Private Member Functions

bool TrackEntersTop (const std::vector< TVector3 > &sorted_trk)
 
bool TrackEntersFront (const std::vector< TVector3 > &sorted_trk)
 
bool TrackEntersBack (const std::vector< TVector3 > &sorted_trk)
 
bool TrackEntersAnode (const std::vector< TVector3 > &sorted_trk, const int driftDir)
 
bool TrackEntersSide (const std::vector< TVector3 > &sorted_trk)
 
bool TrackExitsBottom (const std::vector< TVector3 > &sorted_trk)
 
bool TrackExitsFront (const std::vector< TVector3 > &sorted_trk)
 
bool TrackExitsBack (const std::vector< TVector3 > &sorted_trk)
 
bool TrackExitsAnode (const std::vector< TVector3 > &sorted_trk, const int driftDir)
 
bool TrackExitsSide (const std::vector< TVector3 > &sorted_trk)
 
void SortTrackPoints (const recob::Track &track, std::vector< TVector3 > &sorted_trk)
 
void SplitTrack (const recob::Track &track, std::vector< TVector3 > &sorted_trk)
 
double GetEnteringTimeCoord (const std::vector< TVector3 > &sorted_trk)
 
double GetExitingTimeCoord (const std::vector< TVector3 > &sorted_trk)
 
std::pair< double, size_t > FlashMatch (const double reco_time)
 
double MatchTracks (std::vector< TVector3 > &sorted_trk, std::vector< TLorentzVector > mcpart)
 
std::vector< std::vector< TLorentzVector > > BuildMCParticleList (const art::Handle< std::vector< simb::MCParticle >> &mcpart_h, const double &fTPCResolution, const double &width)
 

Private Attributes

std::string fTrackProducer
 
std::string fT0Producer
 
std::string fFlashProducer
 
std::string fHitProducer
 
std::string fTriggerProducer
 
std::string fCaloProducer
 
bool fUseMC
 
double fTPCResolution
 
double fDriftVelocity
 
bool _debug
 
bool fCathode
 
bool fData
 
double TOP
 
double BOTTOM
 
double FRONT
 
double BACK
 
double det_width
 
std::vector< double > flash_times
 
std::vector< size_t > flash_idx_v
 
double fTimeRes
 
double fPEmin
 
double fRecoT0TimeOffset
 
TTree * tree
 
double mc_time
 
double rc_time
 
double t_match
 
double pe_flash
 
double dt_flash
 
double dt_mc
 
double evttime
 
int year_month_day
 
int hour_min_sec
 
double length
 
double driftDir
 
bool TPC_edge
 
double rc_xs
 
double rc_xe
 
double rc_xs_corr
 
double rc_xe_corr
 
double rc_ys
 
double rc_ye
 
double rc_zs
 
double rc_ze
 
int anode
 
int cathode
 
int sister_track
 
int run
 
int subrun
 
int event
 
int trackNum
 
TTree * evTree
 
int trk_ctr
 
int ev_ctr
 
TTree * caloTree
 
int plane
 
double x
 
double y
 
double z
 
double dQdx
 
double dEdx
 
double dx
 
double resRange
 
double trkLength
 
double kinE
 
double vt0
 

Additional Inherited Members

- Public Types inherited from art::EDAnalyzer
using WorkerType = WorkerT< EDAnalyzer >
 
using ModuleType = EDAnalyzer
 
- Protected Member Functions inherited from art::Observer
std::string const & processName () const
 
bool wantAllEvents () const noexcept
 
bool wantEvent (ScheduleID id, Event const &e) const
 
Handle< TriggerResultsgetTriggerResults (Event const &e) const
 
 Observer (fhicl::ParameterSet const &config)
 
 Observer (std::vector< std::string > const &select_paths, std::vector< std::string > const &reject_paths, fhicl::ParameterSet const &config)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollectorconsumesCollector ()
 
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 50 of file T0RecoSCECalibrations_module.cc.

Constructor & Destructor Documentation

T0RecoSCECalibrations::T0RecoSCECalibrations ( fhicl::ParameterSet const &  p)
explicit

Definition at line 160 of file T0RecoSCECalibrations_module.cc.

161  :
162  EDAnalyzer(p)
163 {
164 
165  fTrackProducer = p.get<std::string>("TrackProducer" );
166  fHitProducer = p.get<std::string>("HitProducer" );
167  fFlashProducer = p.get<std::string>("FlashProducer" );
168  fT0Producer = p.get<std::string>("T0Producer" );
169  fTriggerProducer = p.get<std::string>("TriggerProducer" );
170  fCaloProducer = p.get<std::string>("CaloProducer" );
171  fUseMC = p.get<bool> ("UseMC" );
172  fTPCResolution = p.get<double> ("Resolution" );
173  fTimeRes = p.get<double> ("TimeRes" );
174  fRecoT0TimeOffset = p.get<double> ("RecoT0TimeOffset" );
175  fPEmin = p.get<double> ("PEmin" );
176  fCathode = p.get<bool> ("CathodeOnly" );
177  _debug = p.get<bool> ("debug" );
178  fData = p.get<bool> ("Data" );
179 
180 
181  // get boundaries based on detector bounds
182  auto const* geom = lar::providerFrom<geo::Geometry>();
183 
188 
189  for (geo::TPCID const& tID: geom->IterateTPCIDs()) {
190  geo::TPCGeo const& TPC = geom->TPC(tID);
191 
192  if(TPC.DriftDistance() < 25.0) continue;
193 
194  double origin[3] = {0.};
195  double center[3] = {0.};
196  TPC.LocalToWorld(origin, center);
197 
198  double top = center[1] + TPC.HalfHeight() - fTPCResolution;
199  double bottom = center[1] - TPC.HalfHeight() + fTPCResolution;
200  double front = center[2] - TPC.HalfLength() + fTPCResolution;
201  double back = center[2] + TPC.HalfLength() - fTPCResolution;
202 
203  if (top > TOP) TOP = top;
204  if (bottom < BOTTOM) BOTTOM = bottom;
205  if (front < FRONT) FRONT = front;
206  if (back > BACK) BACK = back;
207 
208  det_width = TPC.DriftDistance();
209  }
210 
211  // Use '_detp' to find 'efield' and 'temp'
212  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();
213  double efield = detProp.Efield();
214  double temp = detProp.Temperature();
215  // Determine the drift velocity from 'efield' and 'temp'
216  fDriftVelocity = detProp.DriftVelocity(efield,temp);
217 }
std::string string
Definition: nybbler.cc:12
Geometry information for a single TPC.
Definition: TPCGeo.h:38
double HalfLength() const
Length is associated with z coordinate [cm].
Definition: TPCGeo.h:117
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
p
Definition: test.py:223
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
double HalfHeight() const
Height is associated with y coordinate [cm].
Definition: TPCGeo.h:111
double DriftDistance() const
Definition: TPCGeo.h:155
def center(depos, point)
Definition: depos.py:117
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:563
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
T0RecoSCECalibrations::T0RecoSCECalibrations ( T0RecoSCECalibrations const &  )
delete
T0RecoSCECalibrations::T0RecoSCECalibrations ( T0RecoSCECalibrations &&  )
delete

Member Function Documentation

void T0RecoSCECalibrations::analyze ( art::Event const &  e)
overridevirtual

Implements art::EDAnalyzer.

Definition at line 298 of file T0RecoSCECalibrations_module.cc.

298  {
299 
300 
301  event = e.event();
302  subrun = e.subRun();
303  run = e.run();
304  trackNum = 0;
305  ev_ctr++;
306 
307  std::vector<std::vector<TLorentzVector>> mcpart_list;
308 
309  if(_debug) std::cout << "NEW EVENT" << std::endl;
310  if (_debug) std::cout << "top: " << TOP << "\nbottom: " << BOTTOM << "\nfront: " << FRONT << "\nback: " << BACK << std::endl;
311 
312 
313  flash_times.clear();
314  flash_idx_v.clear();
315 
316  // load Flash
317  if (_debug) { std::cout << "loading flash from producer " << fFlashProducer << std::endl; }
319  if(!fCathode){
320  flash_h = e.getHandle<std::vector<recob::OpFlash> >(fFlashProducer);
321 
322  // make sure flash looks good
323  if(!flash_h.isValid()) {
324  std::cerr<<"\033[93m[ERROR]\033[00m ... could not locate Flash!"<<std::endl;
325  throw std::exception();
326  }
327  }
328 
329  if(_debug&&fData) std::cout << "loading trigger time from producer " << fTriggerProducer << std::endl;
331  double trigger_time = 0;
332 
333  if(fData){
334  trigger_h = e.getHandle<std::vector<recob::OpFlash> >(fTriggerProducer);
335  trigger_time = trigger_h->at(0).Time();
336  }
337 
338  // load tracks previously created for which T0 reconstruction should occur
339  if (_debug) { std::cout << "loading track from producer " << fTrackProducer << std::endl; }
340  auto track_h = e.getHandle<std::vector<recob::Track> >(fTrackProducer);
341 
342  // make sure tracks look good
343  if(!track_h.isValid()) {
344  std::cerr<<"\033[93m[ERROR]\033[00m ... could not locate Track!"<<std::endl;
345  throw std::exception();
346  }
347 
348  std::vector<art::Ptr<recob::Track> > TrkVec;
349  art::fill_ptr_vector(TrkVec, track_h);
350 
351  // grab 2d hits associated with tracks
352  art::FindMany<recob::Hit> trk_hit_assn_v(track_h, e, fHitProducer);
353 
354  //load calorimetry data products associated with tracks
355  art::FindMany<anab::Calorimetry> trk_calo_assn_v(track_h,e,fCaloProducer);
356 
357  if(_debug){
358  std::cout << "Number of tracks: " << TrkVec.size() << std::endl;
359  std::cout << "Number of hits: " << trk_hit_assn_v.size() << std::endl;
360  std::cout << "Number of calorimetry products: " << trk_calo_assn_v.size() << "\n" << std::endl;
361  }
362 
363  // load MCParticles
364  auto mcpart_h = e.getHandle<std::vector<simb::MCParticle> >("largeant");
365 
366  // if we should use MCParticle
367  if (fUseMC){
368  // make sure particles exist
369  if(!mcpart_h.isValid()) {
370  std::cerr<<"\033[93m[ERROR]\033[00m ... could not locate MCParticle!"<<std::endl;
371  throw std::exception();
372  }
373 
374  mcpart_list = BuildMCParticleList(mcpart_h,fTPCResolution,det_width);
375  }// if use MCParticle
376 
377  // grab T0 objects associated with tracks
378  art::FindMany<anab::T0> trk_t0_assn_v(track_h, e, fT0Producer );
379 
380  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(e);
381  // prepare a vector of optical flash times, if flash above some PE cut value
382  if(!fCathode){
383  size_t flash_ctr = 0;
384  for (auto const& flash : *flash_h){
385  if (flash.TotalPE() > fPEmin){
386  flash_times.push_back( flash.Time() - trigger_time);
387  flash_idx_v.push_back(flash_ctr);
388  if (_debug) std::cout << "\t flash time : " << flash.Time() - trigger_time << ", PE : " << flash.TotalPE() << std::endl;
389  }
390  flash_ctr += 1;
391  }// for all flashes
392 
393  if (_debug) { std::cout << "Selected a total of " << flash_times.size() << " OpFlashes" << std::endl; }
394  }
395 
396  // loop through reconstructed tracks
397  size_t trk_ctr2 = -1;
398 
399  for (auto& track : TrkVec){
400  trk_ctr2 ++;
401  trk_ctr ++;
402  if (_debug) std::cout << "Looping through reco track " << trk_ctr2 << " or " << trk_ctr << std::endl;
403 
404  const std::vector<const recob::Hit*>& Hit_v = trk_hit_assn_v.at(trk_ctr2);
405  if(_debug) std::cout << " Hits loaded" << std::endl;
406  const std::vector<const anab::Calorimetry*>& Calo_v = trk_calo_assn_v.at(trk_ctr2);
407  if(_debug) std::cout << " Calorimetry loaded" << std::endl;
408 
409  mc_time = 0.;
410  rc_time = 0.;
411  t_match = 0.;
412  pe_flash = 0.;
413  dt_flash = -9999.;
414  dt_mc = -9999.;
415  length = 0.;
416  driftDir = 0;
417  rc_xs = 0.;
418  rc_xe = 0.;
419  rc_xs_corr = 0.;
420  rc_xe_corr = 0.;
421  rc_ys = 0.;
422  rc_ye =0.;
423  rc_zs = 0.;
424  rc_ze = 0.;
425  anode = 0;
426  cathode = 0;
427  sister_track = 0;
428  TPC_edge = false;
429 
430  plane = -9999;
431  trkLength = 0.;
432  kinE = 0.;
433  x = 0.;
434  y = 0.;
435  z = 0.;
436  dx = 0.;
437  dEdx = 0.;
438  dQdx = 0.;
439  resRange = 0.;
440  vt0 = 0.;
441 
442  evttime = -9999.;
443  year_month_day = -9999;
444  hour_min_sec = -9999;
445 
446  // get sorted points for the track object [assuming downwards going]
447  std::vector<TVector3> sorted_trk;
448  SortTrackPoints(*track,sorted_trk);
449  if(_debug) std::cout << "\tTrack goes from (" << sorted_trk.at(0).X() << ", " << sorted_trk.at(0).Y() << ", " << sorted_trk.at(0).Z() << ") --> (" << sorted_trk.at(sorted_trk.size()-1).X() << ", " << sorted_trk.at(sorted_trk.size()-1).Y() << ", " << sorted_trk.at(sorted_trk.size()-1).Z() << ")" << std::endl;
450 
451  if( sqrt(pow(sorted_trk.at(0).X() - sorted_trk.at(sorted_trk.size()-1).X(),2.0) + pow(sorted_trk.at(0).Y() - sorted_trk.at(sorted_trk.size()-1).Y(),2.0) + pow(sorted_trk.at(0).Z() - sorted_trk.at(sorted_trk.size()-1).Z(),2.0)) < 50 ){
452  if(_debug) std::cout << "\tTrack too short. Skipping." << std::endl;
453  continue;
454  }
455 
456  // Determine if the track crosses the cathode
457  auto const* geom = lar::providerFrom<geo::Geometry>();
458  auto const* hit = Hit_v.at(0);
459  const geo::WireID wireID = hit->WireID();
460  const auto TPCGeoObject = geom->TPC(wireID.TPC,wireID.Cryostat);
461  short int driftDir1 = TPCGeoObject.DetectDriftDirection();
462  bool cross_cathode = false;
463  for (size_t ii = 1; ii < Hit_v.size(); ii++) {
464  const geo::WireID wireID2 = Hit_v.at(ii)->WireID();
465  const auto TPCGeoObject2 = geom->TPC(wireID2.TPC,wireID2.Cryostat);
466  short int driftDir2 = TPCGeoObject2.DetectDriftDirection();
467 
468  if(driftDir1 + driftDir2 == 0){
469  cross_cathode = true;
470  if(_debug) std::cout << "\tCrosses cathode!" << std::endl;
471  //continue;
472  ii = Hit_v.size();
473  }
474  }
475  if(_debug) std::cout << "\tCross cathode = " << cross_cathode << std::endl;
476 
477 
478  // -------------------------------------------------------------------------------
479  //CATHODE CROSSERS
480  if(cross_cathode){
481  anode = 0;
482  cathode = 1;
483 
484  //const std::vector<const anab::T0*>& T0_v = trk_t0_assn_v.at(trk_ctr);
485  //auto t0 = T0_v.at(0);
486  //rc_time = t0->Time();
487 
488 
489  //Get Time and date
490  art::Timestamp ts = e.time();
491  TTimeStamp tts(ts.timeHigh(), ts.timeLow());
492  evttime = tts.AsDouble();
493 
494  UInt_t year = 0;
495  UInt_t month = 0;
496  UInt_t day = 0;
497  year_month_day = tts.GetDate(kTRUE,0,&year,&month,&day);
498 
499  UInt_t hour = 0;
500  UInt_t min = 0;
501  UInt_t sec = 0;
502  hour_min_sec = tts.GetTime(kTRUE,0,&hour,&min,&sec);
503 
504  SplitTrack(*track,sorted_trk);
505  std::vector<TVector3> top_trk = {sorted_trk.at(0), sorted_trk.at(1)};
506  std::vector<TVector3> bottom_trk = {sorted_trk.at(2), sorted_trk.at(3)};
507 
508  if(_debug)std::cout << "\tCathode-crossing track from (" << top_trk.at(0).X() << ", " << top_trk.at(0).Y() << ", " << top_trk.at(0).Z() << ") --> (" << top_trk.at(1).X() << ", " << top_trk.at(1).Y() << ", " << top_trk.at(1).Z() << ") --> (" << bottom_trk.at(0).X() << ", " << bottom_trk.at(0).Y() << ", " << bottom_trk.at(0).Z() << ") --> (" << bottom_trk.at(1).X() << ", " << bottom_trk.at(1).Y() << ", " << bottom_trk.at(1).Z() << ")" << std::endl;
509 
510  // Top Track!
511  if(_debug) std::cout << "\tTop track" << std::endl;
512 
513  sister_track = trackNum + 1;
514  auto const &top = top_trk.at(0);
515  auto const &bottom = top_trk.at(1);
516 
517  rc_xs = top.X();
518  rc_xs_corr = rc_xs;
519  rc_ys = top.Y();
520  rc_zs = top.Z();
521  rc_xe = bottom.X();
522  rc_xe_corr = rc_xe;
523  rc_ye = bottom.Y();
524  rc_ze = bottom.Z();
525 
526  length = sqrt(pow(rc_xs - rc_xe,2.0) + pow(rc_ys - rc_ye,2.0) + pow(rc_zs - rc_ze,2.0));
527  if(rc_xs<0.0) driftDir = -1;
528  else driftDir = 1;
529 
530  tree->Fill();
531  trackNum++;
532 
533  // Bottom Track
534  if(_debug) std::cout << "\tBottom track" << std::endl;
535  sister_track = trackNum - 1;
536 
537  auto const &top2 = bottom_trk.at(0);
538  auto const &bottom2 = bottom_trk.at(1);
539 
540  rc_xs = top2.X();
541  rc_xs_corr = rc_xs;
542  rc_ys = top2.Y();
543  rc_zs = top2.Z();
544  rc_xe = bottom2.X();
545  rc_xe_corr = rc_xe;
546  rc_ye = bottom2.Y();
547  rc_ze = bottom2.Z();
548 
549  length = sqrt(pow(rc_xs - rc_xe,2.0) + pow(rc_ys - rc_ye,2.0) + pow(rc_zs - rc_ze,2.0));
550  if(rc_xs<0.0) driftDir = -1;
551  else driftDir = 1;
552 
553  tree->Fill();
554 
555 
556  //Save calorimetry information
557  for (size_t ii = 0; ii < Calo_v.size(); ii++) {
558  plane = ii;
559  const anab::Calorimetry* calo_obj = Calo_v.at(ii);
560  trkLength = calo_obj->Range();
561  kinE = calo_obj->KineticEnergy();
562 
563  for (size_t jj = 0; jj < calo_obj->dEdx().size(); jj++){
564  x = calo_obj->XYZ().at(jj).X();
565  y = calo_obj->XYZ().at(jj).Y();
566  z = calo_obj->XYZ().at(jj).Z();
567 
568  dx = calo_obj->TrkPitchVec().at(jj);
569  dEdx = calo_obj->dEdx().at(jj);
570  dQdx = calo_obj->dQdx().at(jj);
571  resRange = calo_obj->ResidualRange().at(jj);
572 
573  caloTree->Fill();
574 
575  }
576 
577  }
578  trackNum++;
579  // ------------------------------------------------------------------------------------
580  // ANODE PIERCERS
581  } else {
582  if(fCathode) {std::cout << "\tSKIPPING ANODE-PEIRCING TRACKS" << std::endl; continue;}
583 
584  if(_debug) std::cout << "\t\tThis track starts in TPC " << wireID.TPC << " which has a drift direction of " << driftDir1 << std::endl;
585  driftDir = driftDir1;
586 
587  // create root trees variables
588  auto const &top = sorted_trk.at(0);
589  auto const &bottom = sorted_trk.at(1);
590 
591  rc_xs = top.X();
592  rc_ys = top.Y();
593  rc_zs = top.Z();
594  rc_xe = bottom.X();
595  rc_ye = bottom.Y();
596  rc_ze = bottom.Z();
597  length = track->Length();
598 
599  //Get Time and date
600  art::Timestamp ts = e.time();
601  TTimeStamp tts(ts.timeHigh(), ts.timeLow());
602  evttime = tts.AsDouble();
603 
604 
605  UInt_t year = 0;
606  UInt_t month = 0;
607  UInt_t day = 0;
608  year_month_day = tts.GetDate(kTRUE,0,&year,&month,&day);
609 
610  UInt_t hour = 0;
611  UInt_t min = 0;
612  UInt_t sec = 0;
613  hour_min_sec = tts.GetTime(kTRUE,0,&hour,&min,&sec);
614 
615 
616  // keep track of whether it goes thorugh the anode or cathode
617  anode = 0, cathode = 0;
618 
619  // 1st category: tracks which ENTER SIDE
620  if ( TrackEntersSide(sorted_trk) == true ) {
621 
622  if (_debug) std::cout << "\t track enters side" << std::endl;
623 
624  // we are not done. We need to check that the track either: 1) Exits the bottom. 2) exits the front or 3) exits the back of the TPC.
625  bool tagged = false;
626 
627  // tracks that exit the bottom
628  if ( (TrackExitsBottom(sorted_trk) == true) ) {
629  tagged = true;
630  if (_debug) std::cout << "\t track exits bottom" << std::endl;
631  }
632  // tracks that exit the front
633  if ( (TrackExitsFront(sorted_trk) == true) and (TrackEntersFront(sorted_trk) == false)) {
634  tagged = true;
635  if (_debug) std::cout << "\t track exits front" << std::endl;
636  }
637  // tracks that exit the back
638  if ( (TrackExitsBack(sorted_trk) == true) and (TrackEntersBack(sorted_trk) == false) ) {
639  tagged = true;
640  if (_debug) std::cout << "\t track exits back" << std::endl;
641  }
642 
643  // has either of these 3 conditions been met? if no, skip this track
644  if (tagged == false) continue;
645 
646  // figure out if it enters the anode or cathode
647  bool enters_anode = TrackEntersAnode(sorted_trk, driftDir);
648 
649  // get the X coordinate of the point piercing the anode/cathode (upon ENTERING)
650  double trkX = GetEnteringTimeCoord(sorted_trk);
651 
652  // reconstruct track T0 w.r.t. trigger time
653 
654  // The 'trkX' enters on the anode, the side of the TPC with a lower x value than the cathode
655  if (enters_anode){
656  rc_time = (det_width - (double)driftDir * trkX) / fDriftVelocity + fRecoT0TimeOffset;
657  if(_debug) std::cout << "\tTrack enters anode. " << rc_time << " = ( " << det_width << " - " << driftDir*trkX << ") / " << fDriftVelocity << std::endl;
658  anode = 1;
659  }
660 
661  }// if the track enters the side
662 
663  // case in which the track exits the side
664  if (TrackExitsSide(sorted_trk) == true) {
665 
666  if (_debug) std::cout << "\t track exits side" << std::endl;
667 
668  // we are not done. We need to check that the track either: 1) Enters the bottom. 2) enters the front or 3) enters the back of the TPC.
669  bool tagged = false;
670 
671  // track enters the top
672  if ( (TrackEntersTop(sorted_trk) == true) ) {
673  tagged = true;
674  if (_debug) std::cout << "\t track enters the top" << std::endl;
675  }
676 
677  if ( (TrackEntersFront(sorted_trk) == true) and (TrackExitsFront(sorted_trk) == false)) {
678  tagged = true;
679  if (_debug) std::cout << "\t track enters front" << std::endl;
680  }
681 
682  if ( (TrackEntersBack(sorted_trk) == true) and (TrackExitsBack(sorted_trk) == false) ) {
683  tagged = true;
684  if (_debug) std::cout << "\t track enters back" << std::endl;
685  }
686 
687  // has either of these 3 conditions been met? if no, skip this track
688  if (tagged == false) continue;
689  // figure out if it enters the anode or cathode
690  bool exits_anode = TrackExitsAnode(sorted_trk, driftDir);
691 
692  // get the X coordinate of the point piercing the anode/cathode (upon ENTERING)
693  double trkX = GetExitingTimeCoord(sorted_trk);
694 
695  // reconstruct track T0 w.r.t. trigger time
696 
697  // The 'trkX' enters on the anode, the side of the TPC with a lower x value than the cathode
698  if (exits_anode){
699  rc_time = (det_width - (double)driftDir * trkX) / fDriftVelocity + fRecoT0TimeOffset;
700  if(_debug) std::cout << "\tTrack exits anode. " << rc_time << " = ( " << det_width << " - " << driftDir*trkX << ") / " << fDriftVelocity << std::endl;
701  anode = 1;
702  }
703 
704  }// if the track exits the side
705 
706  if(!anode) continue;
707 
708 
709  if (_debug)
710  std::cout << "\t this track has a reconstructed time = " << rc_time << std::endl;
711 
714 
715  // Determine if track hits edge of readout window
716  //bool TPC_edge = false;
717  for (auto& hits : Hit_v){
718  auto peakHit = hits->PeakTime();
719  //if(_debug) std::cout << "\t\tHit time track " << trk_ctr << ": " << peakHit << " in TPC " << hits->WireID().TPC << " plane " << hits->WireID().Plane << " and wire " << hits->WireID().Wire << std::endl;
720  if( (peakHit > detProp.ReadOutWindowSize() - 50.0) || (peakHit < 50.0) ){
721  TPC_edge = true;
722  if(_debug) std::cout << "\t\tHit time out of range (0 - " << detProp.ReadOutWindowSize() << "): " << peakHit << std::endl;
723  continue;
724  }
725  }
726 
727  /*if (TPC_edge) {
728  if(_debug) std::cout << "\tHit time too close to edge" << std::endl;
729  continue;
730  }*/
731 
732  // flash matching
733  auto const& flash_match_result = FlashMatch(rc_time);
734  const art::Ptr<recob::OpFlash> flash_ptr(flash_h, flash_match_result.second );
735  if (_debug)
736  std::cout << "\t matched to flash w/ index " << flash_match_result.second << " w/ PE " << flash_ptr->TotalPE() << " and time " << flash_ptr->Time() - trigger_time << " vs reco time " << rc_time << std::endl;
737 
738 
739  pe_flash = flash_ptr->TotalPE();
740  t_match = flash_ptr->Time() - trigger_time;
742 
743  // if we should use MC info -> continue w/ MC validation
744  if (fUseMC == true){
745  // loop through MCParticles to find the one that matches.
746  mc_time = 0;
747  double displacement = 1.0;
748 
749  for (size_t j=0; j < mcpart_list.size(); j++){
750 
751  std::vector<TLorentzVector> mcpart = mcpart_list.at(j);
752 
753  //if(_debug) std::cout <<"\t\tMCParticle: (" << mcpart.at(0).X() << ", " << mcpart.at(0).Y() << ", " << mcpart.at(0).Z() << ") -> (" << mcpart.at(1).X() << ", " << mcpart.at(1).Y() << ", " << mcpart.at(1).Z() << ")" << std::endl;
754 
755  // try matching to MC
756  double tmp_disp = MatchTracks(sorted_trk, mcpart);
757  if ((tmp_disp > displacement)||isnan(tmp_disp))
758  continue;
759 
760  // matched -> get MCTrack time and reconstructed track reconstructed T0
761  displacement = tmp_disp;
762  mc_time = mcpart.at(0).T() / 1000.;
763 
764  if(_debug) std::cout << "\tParticle number: " << j << " with sigma comparison value of " << displacement << "\n\t\tTrack matched to MC value with t0 of " << mc_time << ", flash t0 of " << t_match << " and reconstructed time of " << rc_time << std::endl;
765 
766  } // for all MCParticles
767  }// if we should use MCParticles
768 
769  tree->Fill();
770 
771  //Save calorimetry information
772  for (size_t ii = 0; ii < Calo_v.size(); ii++) {
773  plane = ii;
774  const anab::Calorimetry* calo_obj = Calo_v.at(ii);
775  trkLength = calo_obj->Range();
776  kinE = calo_obj->KineticEnergy();
777  vt0 = driftDir*rc_time*fDriftVelocity;
778 
779  for (size_t jj = 0; jj < calo_obj->dEdx().size(); jj++){
780  x = calo_obj->XYZ().at(jj).X();
781  y = calo_obj->XYZ().at(jj).Y();
782  z = calo_obj->XYZ().at(jj).Z();
783 
784  dx = calo_obj->TrkPitchVec().at(jj);
785  dEdx = calo_obj->dEdx().at(jj);
786  dQdx = calo_obj->dQdx().at(jj);
787  resRange = calo_obj->ResidualRange().at(jj);
788 
789  caloTree->Fill();
790 
791  }
792 
793  }
794 
795  trackNum++;
796 
797  }
798 
799  }
800 
801  evTree->Fill();
802 }
constexpr std::uint32_t timeLow() const
Definition: Timestamp.h:29
double GetEnteringTimeCoord(const std::vector< TVector3 > &sorted_trk)
void SplitTrack(const recob::Track &track, std::vector< TVector3 > &sorted_trk)
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:34
constexpr T pow(T x)
Definition: pow.h:72
bool TrackExitsBack(const std::vector< TVector3 > &sorted_trk)
const std::vector< Point_t > & XYZ() const
Definition: Calorimetry.h:114
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
double MatchTracks(std::vector< TVector3 > &sorted_trk, std::vector< TLorentzVector > mcpart)
std::pair< double, size_t > FlashMatch(const double reco_time)
const std::vector< float > & ResidualRange() const
Definition: Calorimetry.h:103
bool isValid() const noexcept
Definition: Handle.h:191
timescale_traits< TriggerTimeCategory >::time_point_t trigger_time
A point in time on the trigger time scale.
bool TrackExitsAnode(const std::vector< TVector3 > &sorted_trk, const int driftDir)
const double e
const std::vector< float > & dQdx() const
Definition: Calorimetry.h:102
void SortTrackPoints(const recob::Track &track, std::vector< TVector3 > &sorted_trk)
bool TrackEntersSide(const std::vector< TVector3 > &sorted_trk)
bool TrackExitsFront(const std::vector< TVector3 > &sorted_trk)
bool TrackExitsBottom(const std::vector< TVector3 > &sorted_trk)
bool TrackEntersAnode(const std::vector< TVector3 > &sorted_trk, const int driftDir)
std::vector< std::vector< TLorentzVector > > BuildMCParticleList(const art::Handle< std::vector< simb::MCParticle >> &mcpart_h, const double &fTPCResolution, const double &width)
bool TrackEntersTop(const std::vector< TVector3 > &sorted_trk)
Detector simulation of raw signals on wires.
bool TrackEntersBack(const std::vector< TVector3 > &sorted_trk)
const std::vector< float > & dEdx() const
Definition: Calorimetry.h:101
const std::vector< float > & TrkPitchVec() const
Definition: Calorimetry.h:107
bool TrackEntersFront(const std::vector< TVector3 > &sorted_trk)
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
bool TrackExitsSide(const std::vector< TVector3 > &sorted_trk)
const float & KineticEnergy() const
Definition: Calorimetry.h:105
constexpr WireID()=default
Default constructor: an invalid TPC ID.
double GetExitingTimeCoord(const std::vector< TVector3 > &sorted_trk)
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
const float & Range() const
Definition: Calorimetry.h:106
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
QTextStream & endl(QTextStream &s)
void T0RecoSCECalibrations::beginJob ( )
overridevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 219 of file T0RecoSCECalibrations_module.cc.

219  {
220 
222  tree = tfs->make<TTree>("tree","SCE calibrations variables");
223  tree->Branch("mc_time",&mc_time,"mc_time/D");
224  tree->Branch("rc_time",&rc_time,"rc_time/D");
225  tree->Branch("t_match",&t_match,"t_match/D");
226  tree->Branch("dt_flash",&dt_flash,"dt_flash/D");
227  tree->Branch("dt_mc",&dt_flash,"dt_mc/D");
228  tree->Branch("pe_flash",&pe_flash,"pe_flash/D");
229  tree->Branch("length", &length, "length/D");
230  tree->Branch("TPC_edge",&TPC_edge,"TPC_edge/B");
231  tree->Branch("driftDir",&driftDir,"driftDir/I");
232  //Add branches for time of the track
233  tree->Branch("evttime", &evttime, "evttime/D");
234  tree->Branch("year_month_day", &year_month_day, "year_month_day/I");
235  tree->Branch("hour_min_sec", &hour_min_sec, "hour_min_sec/I");
236  // Add branches for the first and last x, y, and z coordinates of the rc tracks and the mc tracks
237  tree->Branch("rc_xs",&rc_xs,"rc_xs/D");
238  tree->Branch("rc_xs_corr",&rc_xs_corr,"rc_xs/D");
239  tree->Branch("rc_ys",&rc_ys,"rc_ys/D");
240  tree->Branch("rc_zs",&rc_zs,"rc_zs/D");
241  tree->Branch("rc_xe",&rc_xe,"rc_xe/D");
242  tree->Branch("rc_xe_corr",&rc_xe_corr,"rc_xe_corr/D");
243  tree->Branch("rc_ye",&rc_ye,"rc_ye/D");
244  tree->Branch("rc_ze",&rc_ze,"rc_ze/D");
245  // information on whether track enters/exits which sides
246  tree->Branch("anode",&anode,"anode/I" );
247  tree->Branch("cathode",&cathode,"cathode/I");
248  tree->Branch("sister_track",&sister_track,"sister_track/I");
249  tree->Branch("run",&run,"run/I");
250  tree->Branch("subrun",&subrun,"subrun/I");
251  tree->Branch("event",&event,"event/I");
252  tree->Branch("trackNum",&trackNum,"trackNum/I");
253 
254  evTree = tfs->make<TTree>("evTree","Event information");
255  evTree->Branch("trk_ctr",&trk_ctr,"trk_ctr/I");
256  evTree->Branch("ev_ctr",&ev_ctr,"ev_ctr/I");
257  ev_ctr = 0;
258  trk_ctr = 0;
259 
260  caloTree = tfs->make<TTree>("caloTree","Uncorrected calorimetry information");
261  caloTree->Branch("anode",&anode,"anode/I" );
262  caloTree->Branch("cathode",&cathode,"cathode/I");
263  caloTree->Branch("run",&run,"run/I");
264  caloTree->Branch("subrun",&subrun,"subrun/I");
265  caloTree->Branch("event",&event,"event/I");
266  caloTree->Branch("trackNum",&trackNum,"trackNum/I");
267  caloTree->Branch("plane",&plane,"plane/I");
268  caloTree->Branch("x",&x,"x/D");
269  caloTree->Branch("y",&y,"y/D");
270  caloTree->Branch("z",&z,"z/D");
271  caloTree->Branch("vt0",&vt0,"vt0/D");
272  caloTree->Branch("dQdx",&dQdx,"dQdx/D");
273  caloTree->Branch("dEdx",&dEdx,"dEdx/D");
274  caloTree->Branch("dx",&dx,"dx/D");
275  caloTree->Branch("resRange",&resRange,"resRange/D");
276  caloTree->Branch("trkLength",&trkLength,"trkLength/D");
277  caloTree->Branch("kinE",&kinE,"kinE/D");
278  caloTree->Branch("dt_flash",&dt_flash,"dt_flash/D");
279  caloTree->Branch("TPC_edge",&TPC_edge,"TPC_edge/B");
280 
281  /*caloSCETree = tfs->make<TTree>("caloTree","SCE corrected calorimetry information");
282  caloSCETree->Branch("run",&run,"run/I");
283  caloSCETree->Branch("subrun",&subrun,"subrun/I");
284  caloSCETree->Branch("event",&event,"event/I");
285  caloSCETree->Branch("trackNum",&trackNum,"trackNum/I");
286  caloSCETree->Branch("x",&x,"x/D");
287  caloSCETree->Branch("y",&y,"y/D");
288  caloSCETree->Branch("z",&z,"z/D");
289  caloSCETree->Branch("dQdx",&dQdx,"dQdx/D");
290  caloSCETree->Branch("dEdx",&dEdx,"dEdx/D");
291  caloSCETree->Branch("dx",&dx,"dx/D");
292  caloSCETree->Branch("resRange",&resRange,"resRange/D");
293  caloSCETree->Branch("trkLength",&trkLength,"trkLength/D");
294  caloSCETree->Branch("kinE",&kinE,"kinE/D");*/
295 
296 }
Event finding and building.
std::vector< std::vector< TLorentzVector > > T0RecoSCECalibrations::BuildMCParticleList ( const art::Handle< std::vector< simb::MCParticle >> &  mcpart_h,
const double &  fTPCResolution,
const double &  width 
)
private

Definition at line 1212 of file T0RecoSCECalibrations_module.cc.

1212  {
1213 
1214  std::vector<std::vector< TLorentzVector>> mcVec;
1215 
1216  for (size_t j=0; j < mcpart_h->size(); j++){
1217 
1218  auto const& mcpart = mcpart_h->at(j);
1219 
1220  double x0 = mcpart.Position(0).X();
1221  double y0 = mcpart.Position(0).Y();
1222  double z0 = mcpart.Position(0).Z();
1223  double t0 = mcpart.T(0);
1224  double x1 = mcpart.Position(mcpart.NumberTrajectoryPoints()-1).X();
1225  double y1 = mcpart.Position(mcpart.NumberTrajectoryPoints()-1).Y();
1226  double z1 = mcpart.Position(mcpart.NumberTrajectoryPoints()-1).Z();
1227  double t1 = mcpart.T(mcpart.NumberTrajectoryPoints()-1);
1228  std::cout << 3 << std::endl;
1229 
1230  //double xs = x0, ys = y0, zs = z0, ts = t0;
1231  //double xe = x1, ye = y1, ze = z1, te = t1;
1232 
1233  bool top_ok = false, bot_ok = false;
1234  std::cout << mcpart.NumberTrajectoryPoints() << std::endl;
1235  if(mcpart.NumberTrajectoryPoints()<2900){
1236  for(size_t k = 0; k<mcpart.NumberTrajectoryPoints(); k++){
1237  if ((x0 > width) || (x0 < -width) || (y0 > TOP + fTPCResolution) || (y0 < BOTTOM - fTPCResolution) || (z0 < FRONT - fTPCResolution) || (z0 > BACK + fTPCResolution) ){
1238  x0 = mcpart.Position(k+1).X();
1239  y0 = mcpart.Position(k+1).Y();
1240  z0 = mcpart.Position(k+1).Z();
1241  t0 = mcpart.T(k+1);
1242  } else top_ok = true;
1243 
1244  if ((x1 > width) || (x1 < -width) || (y1 > TOP + fTPCResolution) || (y1 < BOTTOM - fTPCResolution) || (z1 < FRONT - fTPCResolution) || (z1 > BACK + fTPCResolution) ){
1245  x1 = mcpart.Position(mcpart.NumberTrajectoryPoints()-2-k).X();
1246  y1 = mcpart.Position(mcpart.NumberTrajectoryPoints()-2-k).Y();
1247  z1 = mcpart.Position(mcpart.NumberTrajectoryPoints()-2-k).Z();
1248  t1 = mcpart.T(mcpart.NumberTrajectoryPoints()-2-k);
1249  } else bot_ok = true;
1250 
1251  if(top_ok&&bot_ok) k = mcpart.NumberTrajectoryPoints();
1252  }
1253  }
1254 
1255  double xs = x0, ys = y0, zs = z0, ts = t0;
1256  double xe = x1, ye = y1, ze = z1, te = t1;
1257 
1258  TLorentzVector mc_top(TVector3(xs,ys,zs),ts);
1259  TLorentzVector mc_bot(TVector3(xe,ye,ze),te);
1260  //if (!top_ok||!bot_ok) continue;
1261 
1262  if (ys > ye) mcVec.push_back({mc_top, mc_bot});
1263  else mcVec.push_back({mc_bot, mc_top});
1264  std::cout << 7 << std::endl;
1265  //if(_debug&&_particle) std::cout << " MCParticle: (" << x0 << ", " << y0 << ", " << z0 << ", " << t0/1000. << ") --> (" << xs << ", " << ys << ", " << zs << ", " << ts/1000.0 << ") --> (" << xe << ", " << ye << ", " << ze << ", " << te/1000.0 << ") --> (" << x1 << ", " << y1 << ", " << z1 << ", " << t1/1000.0 << ")" << std::endl;
1266 
1267  }
1268  return mcVec;
1269 
1270 }
code to link reconstructed objects back to the MC truth information
static constexpr double zs
Definition: Units.h:102
static constexpr double ys
Definition: Units.h:103
QTextStream & endl(QTextStream &s)
std::pair< double, size_t > T0RecoSCECalibrations::FlashMatch ( const double  reco_time)
private

Definition at line 867 of file T0RecoSCECalibrations_module.cc.

867  {
868 
869  // loop through all reco'd flash times and see if one matches
870  // the reco time from the track
871  double dt_min = 8000.; // us
872  size_t idx_min = flash_times.size();
873 
874  for (size_t i=0; i < flash_times.size(); i++){
875  auto const& time = flash_times[i];
876  double dt = fabs(time - reco_time);
877  if (dt < dt_min){
878  dt_min = dt;
879  idx_min = flash_idx_v[i];
880  }
881  }
882 
883  std::pair<double,size_t> ret(dt_min,idx_min);
884  return ret;
885 }
double T0RecoSCECalibrations::GetEnteringTimeCoord ( const std::vector< TVector3 > &  sorted_trk)
private

Definition at line 1192 of file T0RecoSCECalibrations_module.cc.

1193 {
1194 
1195  // get the drift-coordinate value
1196  // associated with the point
1197  // along the track piercing the anode / cathode
1198  // ** WHEN the track enters the anode / cathode
1199  return sorted_trk.at(0).X();
1200 }
double T0RecoSCECalibrations::GetExitingTimeCoord ( const std::vector< TVector3 > &  sorted_trk)
private

Definition at line 1203 of file T0RecoSCECalibrations_module.cc.

1204 {
1205  // get the drift-coordinate value
1206  // associated with the point
1207  // along the track piercing the anode / cathode
1208  // ** WHEN the track exits the anode / cathode
1209  return sorted_trk.at(sorted_trk.size() - 1).X();
1210 }
double T0RecoSCECalibrations::MatchTracks ( std::vector< TVector3 > &  sorted_trk,
std::vector< TLorentzVector >  mcpart 
)
private

Definition at line 804 of file T0RecoSCECalibrations_module.cc.

805 {
806  //assumes both the reco track
807  //and mctrack are downwards going
808 
809  //auto const& mctrk_s = mcpart.at(0);
810  //auto const& mctrk_e = mcpart.at(1);
811  //auto const& track_s = sorted_trk.at(0);
812  //auto const& track_e = sorted_trk.at(sorted_trk.size()-1);
813 
814  double mcX = mcpart.at(0).X(), mcY = mcpart.at(0).Y(), mcZ = mcpart.at(0).Z();
815  double mcX_end = mcpart.at(1).X(), mcY_end = mcpart.at(1).Y(), mcZ_end = mcpart.at(1).Z();
816  double trkX = sorted_trk.at(0).X(), trkY = sorted_trk.at(0).Y(), trkZ = sorted_trk.at(0).Z();
817  double trkX_end = sorted_trk.at(sorted_trk.size()-1).X(), trkY_end = sorted_trk.at(sorted_trk.size()-1).Y(), trkZ_end = sorted_trk.at(sorted_trk.size()-1).Z();
818 
819  if(mcY<mcY_end){
820  double tmpX = mcX, tmpY = mcY, tmpZ = mcZ;
821  mcX = mcX_end;
822  mcY = mcY_end;
823  mcZ = mcZ_end;
824  mcX_end = tmpX;
825  mcY_end = tmpY;
826  mcZ_end = tmpZ;
827  }
828 
829  double trkThXZ = atan((trkZ_end - trkZ)/(trkX_end - trkX));
830  double trkThYZ = atan((trkZ_end - trkZ)/(trkY_end - trkY));
831  double sigThXZ = sqrt(2.0)*fTPCResolution/sqrt(pow(trkZ_end-trkZ,2.0)+pow(trkX_end-trkX,2.0));
832  double sigThYZ = sqrt(2.0)*fTPCResolution/sqrt(pow(trkZ_end-trkZ,2.0)+pow(trkY_end-trkY,2.0));
833  double mcThXZ = atan((mcZ_end - mcZ)/(mcX_end - mcX));
834  double mcThYZ = atan((mcZ_end - mcZ)/(mcY_end - mcY));
835 
836  return sqrt(pow((trkThXZ - mcThXZ)/sigThXZ,2.0) + pow((trkThYZ - mcThYZ)/sigThYZ,2.0) + pow((trkY - mcY)/fTPCResolution,2.0) + pow((trkZ - mcZ)/fTPCResolution,2.0))/4.0;
837 
838 
839 
840  /*
841  // if track start is above and mctrk start is above
842  if ( ( track_s.Y() > track_e.Y() ) and ( mctrk_s.Y() > mctrk_e.Y() ) ){
843  if ( (fabs(mctrk_s.Y()-track_s.Y()) < res) and (fabs(mctrk_s.Z()-track_s.Z()) < res) and (fabs(mctrk_e.Y()-track_e.Y()) < res) and (fabs(mctrk_e.Z()-track_e.Z()) < res) )
844  return true;
845  }
846  // if track start is above and mctrk start is below
847  if ( ( track_s.Y() > track_e.Y() ) and ( mctrk_s.Y() < mctrk_e.Y() ) ){
848  if ( (fabs(mctrk_e.Y()-track_s.Y()) < res) and (fabs(mctrk_e.Z()-track_s.Z()) < res) and (fabs(mctrk_s.Y()-track_e.Y()) < res) and (fabs(mctrk_s.Z() - track_e.Z()) < res) )
849  return true;
850  }
851  // if track start is below and mctrk start is above
852  if ( ( track_s.Y() < track_e.Y() ) and ( mctrk_s.Y() > mctrk_e.Y() ) ){
853  if ( (fabs(mctrk_s.Y()-track_e.Y()) < res) and (fabs(mctrk_s.Z()-track_e.Z()) < res) and (fabs(mctrk_e.Y()-track_s.Y()) < res) and (fabs(mctrk_e.Z()-track_s.Z()) < res) )
854  return true;
855  }
856  // if track start is below and mctrk start is below
857  if ( ( track_s.Y() < track_e.Y() ) and ( mctrk_s.Y() < mctrk_e.Y() ) ){
858  if ( (fabs(mctrk_e.Y()-track_e.Y()) < res) and (fabs(mctrk_e.Z()-track_e.Z()) < res) and (fabs(mctrk_s.Y() - track_s.Y()) < res) and (fabs(mctrk_s.Z()-track_s.Z()) < res) )
859  return true;
860  }
861 
862  return false;
863  */
864 }
constexpr T pow(T x)
Definition: pow.h:72
T0RecoSCECalibrations& T0RecoSCECalibrations::operator= ( T0RecoSCECalibrations const &  )
delete
T0RecoSCECalibrations& T0RecoSCECalibrations::operator= ( T0RecoSCECalibrations &&  )
delete
void T0RecoSCECalibrations::SortTrackPoints ( const recob::Track track,
std::vector< TVector3 > &  sorted_trk 
)
private

Definition at line 1136 of file T0RecoSCECalibrations_module.cc.

1137 {
1138  sorted_trk.clear();
1139 
1140  TVector3 track_start, track_end;
1141  double start_y = BOTTOM - 2.0*fTPCResolution;
1142  double end_y = TOP + 2.0*fTPCResolution;
1143 
1144  for (size_t ii = 0; ii < track.NumberTrajectoryPoints(); ii++){
1145  auto const& trk_loc = track.LocationAtPoint(ii);
1146 
1147  if ((trk_loc.X() < -998.)||(trk_loc.Y() < -998.)||(trk_loc.Z() < -998)) continue;
1148 
1149  if (trk_loc.Y() < end_y){
1150  end_y = trk_loc.Y();
1151  track_end = {trk_loc.X(), trk_loc.Y(), trk_loc.Z()};
1152  }
1153  if (trk_loc.Y() > start_y){
1154  start_y = trk_loc.Y();
1155  track_start = {trk_loc.X(), trk_loc.Y(), trk_loc.Z()};
1156  }
1157  }
1158 
1159  sorted_trk.push_back(track_start);
1160  sorted_trk.push_back(track_end);
1161 
1162  /* THIS METHOD ASSUMES THE TRACK IS SORTED AT ALL!!
1163  // vector to store 3D coordinates of
1164  // ordered track
1165  sorted_trk.clear();
1166 
1167  // take the reconstructed 3D track
1168  // and assuming it is downwards
1169  // going, sort points so that
1170  // the track starts at the top
1171  // which point is further up in Y coord?
1172  // start or end?
1173  auto const&N = track.NumberTrajectoryPoints();
1174  auto const&start = track.LocationAtPoint(0);
1175  auto const&end = track.LocationAtPoint( N - 1 );
1176 
1177  // if points are ordered correctly
1178  if (start.Y() > end.Y()){
1179  for (size_t i=0; i < N; i++)
1180  sorted_trk.push_back( track.LocationAtPoint(i) );
1181  }
1182 
1183  // otherwise flip order
1184  else {
1185  for (size_t i=0; i < N; i++)
1186  sorted_trk.push_back( track.LocationAtPoint( N - i - 1) );
1187  }
1188  */
1189 }
Point_t const & LocationAtPoint(size_t i) const
Definition: Track.h:126
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
void T0RecoSCECalibrations::SplitTrack ( const recob::Track track,
std::vector< TVector3 > &  sorted_trk 
)
private

Definition at line 1089 of file T0RecoSCECalibrations_module.cc.

1089  {
1090 
1091  sorted_trk.clear();
1092 
1093  TVector3 track_neg, track_neg_c, track_pos_c,track_pos;
1094  double neg_x = 0.0;
1095  double pos_x = 0.0;
1096  double neg_c = -2.0*fTPCResolution;
1097  double pos_c = 2.0*fTPCResolution;
1098 
1099  for (size_t ii = 0; ii < track.NumberTrajectoryPoints(); ii++){
1100  auto const& trk_loc = track.LocationAtPoint(ii);
1101 
1102  if ((trk_loc.X() < -998.)||(trk_loc.Y() < -998.)||(trk_loc.Z() < -998)) continue;
1103 
1104  if (trk_loc.X() < neg_x){
1105  neg_x = trk_loc.X();
1106  track_neg = {trk_loc.X(), trk_loc.Y(), trk_loc.Z()};
1107  }
1108  if (trk_loc.X() > pos_x){
1109  pos_x = trk_loc.X();
1110  track_pos = {trk_loc.X(), trk_loc.Y(), trk_loc.Z()};
1111  }
1112  if ((trk_loc.X() < 0.0) && (trk_loc.X() > neg_c)){
1113  neg_c = trk_loc.X();
1114  track_neg_c = {trk_loc.X(), trk_loc.Y(), trk_loc.Z()};
1115  }
1116  if ((trk_loc.X() > 0.0) && (trk_loc.X() < pos_c)){
1117  pos_c = trk_loc.X();
1118  track_pos_c = {trk_loc.X(), trk_loc.Y(), trk_loc.Z()};
1119  }
1120  }
1121 
1122  if( track_neg.Y() > track_pos.Y()){
1123  sorted_trk.push_back(track_neg);
1124  sorted_trk.push_back(track_neg_c);
1125  sorted_trk.push_back(track_pos_c);
1126  sorted_trk.push_back(track_pos);
1127  } else {
1128  sorted_trk.push_back(track_pos);
1129  sorted_trk.push_back(track_pos_c);
1130  sorted_trk.push_back(track_neg_c);
1131  sorted_trk.push_back(track_neg);
1132  }
1133 
1134 }
Point_t const & LocationAtPoint(size_t i) const
Definition: Track.h:126
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
bool T0RecoSCECalibrations::TrackEntersAnode ( const std::vector< TVector3 > &  sorted_trk,
const int  driftDir 
)
private

Definition at line 946 of file T0RecoSCECalibrations_module.cc.

947 {
948 
949  // we know the track enters either the
950  // anode or cathode
951  // at this point figure out
952  // if it ENTERS the ANODE or CATHODE
953  // ANODE: top point must be at lower X-coord
954  // than bottom point
955  // CATHODE: top point must be at larger X-coord
956  // than bottom point
957  // assume track has already been sorted
958  // such that the 1st point is the most elevated in Y coord.
959  // return TRUE if passes the ANODE
960 
961  auto const& top = sorted_trk.at(0);
962  auto const& bottom = sorted_trk.at( sorted_trk.size() - 1 );
963 
964  if ( ( (top.X() < bottom.X()) && driftDir<0 ) || ( (top.X() > bottom.X()) && driftDir>0 ) )
965  return true;
966 
967  return false;
968 }
bool T0RecoSCECalibrations::TrackEntersBack ( const std::vector< TVector3 > &  sorted_trk)
private

Definition at line 925 of file T0RecoSCECalibrations_module.cc.

926 {
927 
928  // Determines if the track enters the
929  // back of the TPC based on if the position
930  // of its initial Z-coordinate is greater
931  // than the location of the back of the
932  // TPC in Z
933 
934  // First define 'top_pt' to mean the point at the start of the track
935  auto const& top_pt = sorted_trk.at(0);
936 
937  if (top_pt.Z() > BACK)
938  return true;
939 
940  // If this condition is not satisfied, then return 'false' (the track was not determined
941  // within resolution to enter the back of the TPC)
942  return false;
943 }
bool T0RecoSCECalibrations::TrackEntersFront ( const std::vector< TVector3 > &  sorted_trk)
private

Definition at line 901 of file T0RecoSCECalibrations_module.cc.

902 {
903 
904  // Determine if the track enters the
905  // front of the TPC based on if the position
906  // of its initial Z-coordinate is less than
907  // the location of the front of the TPC in Z
908 
909  // First define 'top_pt' to mean the point at the start of the track
910  auto const& top_pt = sorted_trk.at(0);
911 
912  if (top_pt.Z() < FRONT)
913  return true;
914 
915  // I may include the case in which I check
916  // the y-coordinates as well, but I will not
917  // implement that at this time
918 
919  // If this condition is not satisfied, then return 'false' (the track was not determined
920  // within resolution to enter the front of the TPC)
921  return false;
922 }
bool T0RecoSCECalibrations::TrackEntersSide ( const std::vector< TVector3 > &  sorted_trk)
private

Definition at line 971 of file T0RecoSCECalibrations_module.cc.

972 {
973 
974  // check that the top-most point
975  // is not on the top of the TPC
976  // nor on the front & back of the TPC
977 
978  auto const& top_pt = sorted_trk.at(0);
979 
980  // if highest point above the TOP -> false
981  if (top_pt.Y() > TOP)
982  return false;
983 
984  // if highest point in Z close to front or back
985  // -> FALSE
986  if ( (top_pt.Z() < FRONT) or (top_pt.Z() > BACK) )
987  return false;
988 
989 
990  // If the function makes it this far, then it will enter through one of the sides of the TPC
991  return true;
992 }
bool T0RecoSCECalibrations::TrackEntersTop ( const std::vector< TVector3 > &  sorted_trk)
private

Definition at line 888 of file T0RecoSCECalibrations_module.cc.

889 {
890  // check that the first point in the track
891  // pierces the top boundary of the TPC
892  // This track either will pierce the top of the TPC or is just about to (the '_TOP' variable is just below the actual coordinate position of the top in Y)
893 
894  if (sorted_trk.at(0).Y() > TOP)
895  return true;
896 
897  return false;
898 }
bool T0RecoSCECalibrations::TrackExitsAnode ( const std::vector< TVector3 > &  sorted_trk,
const int  driftDir 
)
private

Definition at line 1043 of file T0RecoSCECalibrations_module.cc.

1044 {
1045 
1046  // Check, once it's known that the track doesn't exit out of the bottom, whether it's the anode or
1047  // the cathode that it exits out of
1048  // This can be done by direct analogy with the 'Anode' function (shown in this file as the 'TrackEntersAnode') function written by D. Caratelli
1049  // Define 'top' as the point at the start of the track, and 'bottom' as the point at the end of the track
1050 
1051  auto const& top = sorted_trk.at(0);
1052  auto const& bottom = sorted_trk.at(sorted_trk.size() - 1);
1053 
1054  // Check to see which point has a lower x coordinate
1055  // If the bottom does, then it exits out of the anode
1056  // If the top does, then it exits out of the cathode
1057  if ( ( (bottom.X() < top.X()) && driftDir<0 ) || ( (bottom.X() > top.X()) && driftDir>0 ) )
1058  return true;
1059 
1060  return false; // Otherwise, the top is less than the bottom, so the track ended closer to the cathode and exited there
1061 }
bool T0RecoSCECalibrations::TrackExitsBack ( const std::vector< TVector3 > &  sorted_trk)
private

Definition at line 1025 of file T0RecoSCECalibrations_module.cc.

1026 {
1027 
1028  // Determine if the track exits the
1029  // front of the TPC based on if the position
1030  // of its final Z-coordinate is less than
1031  // the location of the front of the TPC in Z
1032 
1033  // First define 'bottom_pt' to mean the point at the end of the track
1034  auto const& bottom_pt = sorted_trk.at(sorted_trk.size() - 1);
1035 
1036  if (bottom_pt.Z() > BACK)
1037  return true;
1038 
1039  return false;
1040 }
bool T0RecoSCECalibrations::TrackExitsBottom ( const std::vector< TVector3 > &  sorted_trk)
private

Definition at line 995 of file T0RecoSCECalibrations_module.cc.

996 {
997 
998  // check that the last point in the track
999  // pierces the bottom boundary of the TPC
1000  if ( sorted_trk.at( sorted_trk.size() - 1).Y() < BOTTOM )
1001  return true;
1002 
1003  return false;
1004 }
bool T0RecoSCECalibrations::TrackExitsFront ( const std::vector< TVector3 > &  sorted_trk)
private

Definition at line 1007 of file T0RecoSCECalibrations_module.cc.

1008 {
1009 
1010  // Determine if the track exits the
1011  // front of the TPC based on if the position
1012  // of its final Z-coordinate is less than
1013  // the location of the front of the TPC in Z
1014 
1015  // First define 'bottom_pt' to mean the point at the end of the track
1016  auto const& bottom_pt = sorted_trk.at(sorted_trk.size() - 1);
1017 
1018  if (bottom_pt.Z() < FRONT)
1019  return true;
1020 
1021  return false;
1022 }
bool T0RecoSCECalibrations::TrackExitsSide ( const std::vector< TVector3 > &  sorted_trk)
private

Definition at line 1064 of file T0RecoSCECalibrations_module.cc.

1065 {
1066 
1067  // check that the bottom-most point
1068  // is not on the bottom of the TPC
1069  // nor on the front & back of the TPC
1070 
1071  auto const& bottom_pt = sorted_trk.at(sorted_trk.size() - 1);
1072 
1073  // if lowest point below the BOTTOM -> false
1074  // Within this resolution, this means that it's likely that the track exited out of the bottom (at a point earlier on in the process than the last point) OR is just about to
1075 
1076  if (bottom_pt.Y() < BOTTOM)
1077  return false;
1078 
1079  // if lowest point in Z close to front or back
1080  // -> FALSE
1081  // If the the bottom point is less than the front, then the track has already pierced the front of the TPC and exited that way OR is likely just about to
1082  // If the bottom point is greater than the back, then the track has already pierced the back of the TPC and exited that way OR is likely just about to
1083  if ( (bottom_pt.Z() < FRONT) or (bottom_pt.Z() > BACK) )
1084  return false;
1085 
1086  return true;
1087 }

Member Data Documentation

bool T0RecoSCECalibrations::_debug
private

Definition at line 81 of file T0RecoSCECalibrations_module.cc.

int T0RecoSCECalibrations::anode
private

Definition at line 139 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::BACK
private

Definition at line 85 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::BOTTOM
private

Definition at line 85 of file T0RecoSCECalibrations_module.cc.

TTree* T0RecoSCECalibrations::caloTree
private

Definition at line 149 of file T0RecoSCECalibrations_module.cc.

int T0RecoSCECalibrations::cathode
private

Definition at line 140 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::dEdx
private

Definition at line 153 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::det_width
private

Definition at line 85 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::dQdx
private

Definition at line 153 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::driftDir
private

Definition at line 133 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::dt_flash
private

Definition at line 125 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::dt_mc
private

Definition at line 126 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::dx
private

Definition at line 153 of file T0RecoSCECalibrations_module.cc.

int T0RecoSCECalibrations::ev_ctr
private

Definition at line 147 of file T0RecoSCECalibrations_module.cc.

int T0RecoSCECalibrations::event
private

Definition at line 143 of file T0RecoSCECalibrations_module.cc.

TTree* T0RecoSCECalibrations::evTree
private

Definition at line 145 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::evttime
private

Definition at line 128 of file T0RecoSCECalibrations_module.cc.

std::string T0RecoSCECalibrations::fCaloProducer
private

Definition at line 75 of file T0RecoSCECalibrations_module.cc.

bool T0RecoSCECalibrations::fCathode
private

Definition at line 82 of file T0RecoSCECalibrations_module.cc.

bool T0RecoSCECalibrations::fData
private

Definition at line 83 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::fDriftVelocity
private

Definition at line 79 of file T0RecoSCECalibrations_module.cc.

std::string T0RecoSCECalibrations::fFlashProducer
private

Definition at line 72 of file T0RecoSCECalibrations_module.cc.

std::string T0RecoSCECalibrations::fHitProducer
private

Definition at line 73 of file T0RecoSCECalibrations_module.cc.

std::vector<size_t> T0RecoSCECalibrations::flash_idx_v
private

Definition at line 88 of file T0RecoSCECalibrations_module.cc.

std::vector<double> T0RecoSCECalibrations::flash_times
private

Definition at line 87 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::fPEmin
private

Definition at line 92 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::fRecoT0TimeOffset
private

Definition at line 94 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::FRONT
private

Definition at line 85 of file T0RecoSCECalibrations_module.cc.

std::string T0RecoSCECalibrations::fT0Producer
private

Definition at line 71 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::fTimeRes
private

Definition at line 90 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::fTPCResolution
private

Definition at line 78 of file T0RecoSCECalibrations_module.cc.

std::string T0RecoSCECalibrations::fTrackProducer
private

Definition at line 70 of file T0RecoSCECalibrations_module.cc.

std::string T0RecoSCECalibrations::fTriggerProducer
private

Definition at line 74 of file T0RecoSCECalibrations_module.cc.

bool T0RecoSCECalibrations::fUseMC
private

Definition at line 77 of file T0RecoSCECalibrations_module.cc.

int T0RecoSCECalibrations::hour_min_sec
private

Definition at line 130 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::kinE
private

Definition at line 154 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::length
private

Definition at line 132 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::mc_time
private

Definition at line 121 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::pe_flash
private

Definition at line 124 of file T0RecoSCECalibrations_module.cc.

int T0RecoSCECalibrations::plane
private

Definition at line 151 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::rc_time
private

Definition at line 122 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::rc_xe
private

Definition at line 135 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::rc_xe_corr
private

Definition at line 135 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::rc_xs
private

Definition at line 135 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::rc_xs_corr
private

Definition at line 135 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::rc_ye
private

Definition at line 136 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::rc_ys
private

Definition at line 136 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::rc_ze
private

Definition at line 137 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::rc_zs
private

Definition at line 137 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::resRange
private

Definition at line 153 of file T0RecoSCECalibrations_module.cc.

int T0RecoSCECalibrations::run
private

Definition at line 143 of file T0RecoSCECalibrations_module.cc.

int T0RecoSCECalibrations::sister_track
private

Definition at line 141 of file T0RecoSCECalibrations_module.cc.

int T0RecoSCECalibrations::subrun
private

Definition at line 143 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::t_match
private

Definition at line 123 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::TOP
private

Definition at line 85 of file T0RecoSCECalibrations_module.cc.

bool T0RecoSCECalibrations::TPC_edge
private

Definition at line 134 of file T0RecoSCECalibrations_module.cc.

int T0RecoSCECalibrations::trackNum
private

Definition at line 143 of file T0RecoSCECalibrations_module.cc.

TTree* T0RecoSCECalibrations::tree
private

Definition at line 120 of file T0RecoSCECalibrations_module.cc.

int T0RecoSCECalibrations::trk_ctr
private

Definition at line 146 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::trkLength
private

Definition at line 154 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::vt0
private

Definition at line 155 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::x
private

Definition at line 152 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::y
private

Definition at line 152 of file T0RecoSCECalibrations_module.cc.

int T0RecoSCECalibrations::year_month_day
private

Definition at line 129 of file T0RecoSCECalibrations_module.cc.

double T0RecoSCECalibrations::z
private

Definition at line 152 of file T0RecoSCECalibrations_module.cc.


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