Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions | Private Member Functions | Private Attributes | List of all members
upmuana::UpMuAnalysis Class Reference
Inheritance diagram for upmuana::UpMuAnalysis:
art::EDAnalyzer art::EventObserverBase art::Consumer art::EngineCreator

Public Types

using WorkerType = WorkerT< EDAnalyzer >
 
using ModuleType = EDAnalyzer
 

Public Member Functions

 UpMuAnalysis (fhicl::ParameterSet const &p)
 
 UpMuAnalysis (UpMuAnalysis const &)=delete
 
 UpMuAnalysis (UpMuAnalysis &&)=delete
 
UpMuAnalysisoperator= (UpMuAnalysis const &)=delete
 
UpMuAnalysisoperator= (UpMuAnalysis &&)=delete
 
void analyze (art::Event const &e) override
 
void beginJob () override
 
std::string workerType () const
 
bool modifiesEvent () const
 
void registerProducts (MasterProductRegistry &, ProductDescriptions &, ModuleDescription const &)
 
std::string const & processName () const
 
bool wantAllEvents () const
 
bool wantEvent (Event const &e)
 
fhicl::ParameterSetID selectorConfig () const
 
art::Handle< art::TriggerResultsgetTriggerResults (Event const &e) const
 
template<typename T , BranchType = InEvent>
ProductToken< Tconsumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< Tconsumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TconsumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< TmayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< TmayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TmayConsumeView (InputTag const &it)
 
base_engine_tcreateEngine (seed_t seed)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumernon_module_context ()
 

Protected Member Functions

CurrentProcessingContext const * currentContext () const
 
detail::CachedProducts & cachedProducts ()
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

Private Member Functions

float containmentType (rb::Track)
 
double getLLR (std::set< std::pair< rb::CellHit, rb::RecoHit >, bool(*)(std::pair< rb::CellHit, rb::RecoHit >, std::pair< rb::CellHit, rb::RecoHit >)>, std::vector< rb::RecoHit > &outliers, double &P_up, double &P_dn, double &chi2, double &slope)
 
void LLR (std::vector< double > &eT, std::vector< double > &mT, std::vector< double > &mTErr, double &slope, double &chi2, double &P_up, double &P_dn, std::vector< rb::RecoHit > &in, std::vector< rb::RecoHit > &outliers)
 
void LinFit (const std::vector< double > &x, const std::vector< double > &y, double *fitpar)
 
void LinFit (const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &ye, double *fitpar)
 
TVector3 AnglesToVector (double zen, double azi) const
 

Private Attributes

art::InputTag trackTag_
 
art::InputTag fSliceModuleLabel
 
art::InputTag fSliceInstance
 
double fMinY = -700.0
 
double fMinX = -650.0
 
double fMinZ = 100.0
 
double fMaxY = 400.0
 
double fMaxX = 650.0
 
double fMaxZ = 5800.0
 
art::ServiceHandle< locator::CelestialLocatorfSunPos
 
trk::CosmicTrackUtilities fUtilities
 
TNtuple * ntp_track
 

Detailed Description

Definition at line 77 of file UpMuAnalysis_module.cc.

Member Typedef Documentation

Definition at line 39 of file EDAnalyzer.h.

Definition at line 38 of file EDAnalyzer.h.

Constructor & Destructor Documentation

upmuana::UpMuAnalysis::UpMuAnalysis ( fhicl::ParameterSet const &  p)
explicit

Definition at line 134 of file UpMuAnalysis_module.cc.

135  :
136  EDAnalyzer(p) // ,
137  // More initializers here.
138  , trackTag_(p.get<std::string>("trackInputTag"))
139  , fSliceModuleLabel (p.get<std::string>("SliceModuleLabel"))
140  , fSliceInstance (p.get<std::string>("SliceInstanceLabel"))
141  , ntp_track(nullptr)
142 {}
const char * p
Definition: xmltok.h:285
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
enum BeamMode string
upmuana::UpMuAnalysis::UpMuAnalysis ( UpMuAnalysis const &  )
delete
upmuana::UpMuAnalysis::UpMuAnalysis ( UpMuAnalysis &&  )
delete

Member Function Documentation

void upmuana::UpMuAnalysis::analyze ( art::Event const &  e)
overridevirtual

Implements art::EDAnalyzer.

Definition at line 152 of file UpMuAnalysis_module.cc.

References std::acos(), rb::Cluster::AllCells(), AnglesToVector(), art::PtrVector< T >::at(), std::atan(), std::atan2(), trk::CosmicTrackUtilities::Azimuth(), rb::CellHit::Cell(), chi2(), containmentType(), novadaq::timeutils::convertNovaTimeToUnixTime(), trk::CosmicTrackUtilities::CosTheta(), art::errors::DataCorruption, rb::Prong::Dir(), dist, febshutoff_auto::end, allTimeWatchdog::endl, art::EventID::event(), check_time_usage::float, fSunPos, fUtilities, art::DataViewImpl::getByLabel(), getLLR(), locator::CelestialLocator::GetSunPosition_FD(), rb::RecoHit::GeV(), rb::CellHit::GoodTiming(), h1, h2, MECModelEnuComparisons::i, art::Event::id(), makeTrainCVSamples::int, rb::Track::InterpolateDir(), rb::Prong::Is3D(), rb::RecoHit::IsCalibrated(), geo::kX, geo::kY, util::LinFit(), LinFit(), M_PI, make_pair(), ntp_track, upmuana::pairHitComp(), rb::CellHit::Plane(), art::PtrVector< T >::push_back(), rb::Cluster::RecoHit(), art::EventID::run(), moon_position_table_new3::second, art::PtrVector< T >::size(), std::sqrt(), febshutoff_auto::start, rb::Prong::Start(), rb::Track::Stop(), art::EventID::subRun(), std::swap(), rb::RecoHit::T(), msf_helper::timespec, rb::Track::TotalLength(), rawdata::DAQHeader::TotalMicroSlices(), Unit(), POTSpillRate::view, rb::CellHit::View(), wgt, rb::RecoHit::X(), rb::Cluster::XCells(), rb::RecoHit::Y(), rb::Cluster::YCells(), and rb::RecoHit::Z().

153 {
154  // Get all slices
156  e.getByLabel( fSliceModuleLabel, slicecol);
158 
159  for(int i = 0; i < (int)slicecol->size(); i++){
160  art::Ptr<rb::Cluster> sli(slicecol, i);
161  slces.push_back(sli);
162  }
163 
164  // get all
166  e.getByLabel(trackTag_, tracks);
167  art::FindOneP<rb::Cluster> fmSlice(tracks, e, trackTag_);
168 
170  e.getByLabel("daq", raw_triggers); // this is need it
171  if (raw_triggers->size() != 1)
173  << __FILE__ << ":" << __LINE__ << "\n"
174  << "RawTrigger vector has incorrect size: " << raw_triggers->size()
175  << std::endl;
176 
177  std::vector<const rb::Cluster*> slices;
178  std::vector<float> slice_containment;
179 
180  art::Handle<rawdata::DAQHeader> raw_daqheader;
181  e.getByLabel("daq", raw_daqheader);
182 
183  double zen, azi;
184  unsigned long long event_time(raw_triggers->front().fTriggerTimingMarker_TimeStart);
185 
186  struct timespec ts;
188  convertNovaTimeToUnixTime(event_time, ts);
189 
190  fSunPos->GetSunPosition_FD(ts.tv_sec, zen, azi);
191 
192  //double zenR, aziR;
193  //fSunPos->GetRndmSunPosition_FD(ts, zenR, aziR);
194 
195  std::vector<unsigned> tracks_per_slice;
196  for (size_t i_track=0; i_track < tracks->size(); ++i_track){
197 
198  rb::Track theTrack = tracks->at(i_track);
199  float containment = containmentType(theTrack);
200 
201  float slice=-1;
202  bool found_slice = false;
203 
204  // Loop over Slices
205  const rb::Cluster *theSlice = &(*fmSlice.at(i_track));
206 
207  for (size_t i_slice=0; i_slice<slices.size(); ++i_slice) {
208  if (theSlice == slices.at(i_slice)){
209  slice = (float)i_slice; found_slice = true;
210  ++tracks_per_slice.at(i_slice);
211  //do the linear fit on the slice
212  if(containment != 4) slice_containment.at(i_slice) = 1;
213  break;
214  }
215  } //end loop slice
216 
217  if (!found_slice){
218  slice = (float)slices.size();
219  slices.push_back(theSlice);
220  tracks_per_slice.push_back(1);
221  if(containment!=4)
222  slice_containment.push_back(1);
223  else
224  slice_containment.push_back(4);
225  }
226 
227  // another loop to find slice id from vector slices
228  slice = -1;
229  for(size_t i_slice = 0; i_slice < slces.size(); i_slice++)
230  {
231  if(slces.at(i_slice).get() == theSlice)
232  {
233  slice = (float)i_slice;
234  break;
235  }
236  }
237 
238 
239  if (!theTrack.Is3D()) continue;
240 
241  art::PtrVector<rb::CellHit> trackHits = theTrack.AllCells();
242  std::set< std::pair<rb::CellHit,rb::RecoHit>,
243  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
244  std::pair<rb::CellHit,rb::RecoHit>)>
245  sortedTrackHits (pairHitComp);
246  std::set< std::pair<rb::CellHit,rb::RecoHit>,
247  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
248  std::pair<rb::CellHit,rb::RecoHit>)>
249  sortedTrackHitsX (pairHitComp);
250  std::set< std::pair<rb::CellHit,rb::RecoHit>,
251  bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
252  std::pair<rb::CellHit,rb::RecoHit>)>
253  sortedTrackHitsY (pairHitComp);
254 
255  std::vector<double> x_hit;
256  std::vector<double> y_hit;
257  std::vector<double> zy_hit;
258  std::vector<double> zx_hit;
259 
260  // Loop over Hits
261  for (size_t i_hit=0; i_hit < trackHits.size(); ++i_hit) {
262  art::Ptr<rb::CellHit> theHit = trackHits.at(i_hit);
263  if (!theHit->GoodTiming()) continue;
264  rb::RecoHit theRecoHit = theTrack.RecoHit(theHit);
265  if (!theRecoHit.IsCalibrated()) continue;
266  sortedTrackHits.insert(std::make_pair(*theHit,theRecoHit));
267 
268  unsigned short iC = theHit->Cell();
269  unsigned short iP = theHit->Plane();
270  geo::View_t view = theHit->View();
271  if(view == geo::kY){
272  x_hit.push_back(iC);
273  zx_hit.push_back(iP);
274  sortedTrackHitsY.insert(std::make_pair(*theHit,theRecoHit));
275  }
276  else if(view == geo::kX){
277  y_hit.push_back(iC);
278  zy_hit.push_back(iP);
279  sortedTrackHitsX.insert(std::make_pair(*theHit,theRecoHit));
280  }
281  } // end loop hit
282 
283  double fitpar[3];
284  LinFit(x_hit, zx_hit, fitpar);
285  double r2x = fitpar[2];
286  LinFit(y_hit, zy_hit, fitpar);
287  double r2y = fitpar[2];
288  unsigned nxhit = x_hit.size();
289  unsigned nyhit = y_hit.size();
290 
291  //////////////////////////////////////////////////////////////////////////
292  //cp linear fit
293 
294  Double_t chi2X_fit = 1e6;
295  Double_t chi2Y_fit = 1e6;
296  const double wgt = 1;
297  Double_t trackX_fromfit = 0;
298  Double_t trackY_fromfit = 0;
299 
300  Double_t trackX_fromfit_slice = 0;
301  Double_t trackY_fromfit_slice = 0;
302 
303  double mvalX, cvalX, mvalY, cvalY;
304  TH1F *h1 = new TH1F("h1","xline",50,0,600);
305  TH1F *h2 = new TH1F("h2","yline",50,0,600);
306  // Loop over Slices
307  //for (size_t i = 0; i<slices.size(); ++i) {
308  std::vector<double> x_hit_fit;
309  std::vector<double> y_hit_fit;
310  std::vector<double> zy_hit_fit;
311  std::vector<double> zx_hit_fit;
312 
313  std::vector<double> weight_x;
314  std::vector<double> weight_y;
315 
316 
317  //const rb::Cluster* slice = theSlice;
318 
319  art::PtrVector<rb::CellHit> xcells = theSlice->XCells();
320  art::PtrVector<rb::CellHit> ycells = theSlice->YCells();
321 
322  for(size_t y_ind =0; y_ind < ycells.size(); y_ind++)
323  {
324  art::Ptr<rb::CellHit> theHity = ycells.at(y_ind);
325 
326  unsigned short iC = theHity->Cell();
327  unsigned short iP = theHity->Plane();
328  y_hit_fit.push_back(iC);
329  zy_hit_fit.push_back(iP);
330  weight_y.push_back(wgt);
331  trackY_fromfit += y_ind;
332  h1->Fill(y_hit_fit.at(y_ind));
333  }
334 
335  for(size_t x_ind = 0; x_ind < xcells.size(); x_ind++)
336  {
337  art::Ptr<rb::CellHit> theHitx = xcells.at(x_ind);
338 
339  unsigned short iC = theHitx->Cell();
340  unsigned short iP = theHitx->Plane();
341  x_hit_fit.push_back(iC);
342  zx_hit_fit.push_back(iP);
343  weight_x.push_back(wgt);
344  trackX_fromfit += x_ind;
345  h2->Fill(x_hit_fit.at(x_ind));
346  }
347  trackY_fromfit_slice += trackY_fromfit;
348  trackX_fromfit_slice += trackX_fromfit;
349  //}
350  if(x_hit_fit.size()>=2){
351  double chi2X_aux = util::LinFit(zx_hit_fit, x_hit_fit, weight_x, mvalX, cvalX);
352  if(chi2X_aux==chi2X_aux)chi2X_fit = chi2X_aux;
353  }
354  if(y_hit_fit.size()>=2){
355  double chi2Y_aux = util::LinFit(zy_hit_fit, y_hit_fit, weight_y, mvalY, cvalY);
356  if(chi2Y_aux==chi2Y_aux)chi2Y_fit = chi2Y_aux;
357  }
358 
359  ///////////////////////////////////////////////////////////
360  std::vector<rb::RecoHit> outliers, outliersX, outliersY;
361  double llr, P_up, P_dn, chi2, slope;
362  double llrX, P_upX, P_dnX, chi2X, slopeX;
363  double llrY, P_upY, P_dnY, chi2Y, slopeY;
364 
365  llr = getLLR(sortedTrackHits, outliers, P_up, P_dn, chi2, slope);
366  llrX = getLLR(sortedTrackHitsX, outliersX, P_upX, P_dnX, chi2X, slopeX);
367  llrY = getLLR(sortedTrackHitsY, outliersY, P_upY, P_dnY, chi2Y, slopeY);
368 
371  if (sortedTrackHits.size() >= 1) start = (sortedTrackHits.begin()->second);
372  else { start = rb::RecoHit(); end = rb::RecoHit(); }
373  if (sortedTrackHits.size() >= 2) end = (--sortedTrackHits.end())->second;
374  else end = start;
375  float dist = theTrack.TotalLength();
376 
377  // Average values
378  float tDirX=0, tDirY=0, tDirZ=0, tEleAngle, trackTotalE=0;
379  float avgTX=0, avgTY=0, tCosTheta, tAzimuth;
380 
381  for (size_t i_hit=0; i_hit < trackHits.size(); ++i_hit) {
382  art::Ptr<rb::CellHit> theHit = trackHits.at(i_hit);
383  rb::RecoHit theRecoHit = theTrack.RecoHit(theHit);
384 
385  TVector3 theDir = theTrack.InterpolateDir(theRecoHit.Z());
386 
387  tDirX += theDir.x()/trackHits.size();
388  tDirY += theDir.y()/trackHits.size();
389  tDirZ += theDir.z()/trackHits.size();
390 
391  if(theRecoHit.IsCalibrated()) trackTotalE += theRecoHit.GeV();
392 
393  if (theHit->View() == geo::kX)
394  avgTX += theRecoHit.T()/sortedTrackHitsX.size();
395  else if (theHit->View() == geo::kY)
396  avgTY += theRecoHit.T()/sortedTrackHitsY.size();
397  } // close the Hit loop
398 
399  tEleAngle = atan(tDirY/sqrt(tDirX*tDirX + tDirZ*tDirZ));
400  tCosTheta = fUtilities.CosTheta(theTrack.Dir().Y(), theTrack.Dir().Mag());
401  tAzimuth = fUtilities.Azimuth(theTrack.Dir().X(), theTrack.Dir().Z(), theTrack.Dir().Mag());
402 
403  TVector3 start_trk = theTrack.Start();
404  TVector3 end_trk = theTrack.Stop();
405  if(start_trk.Y() < end_trk.Y()) std::swap(start_trk, end_trk);
406  const TVector3 dir_trk = (start_trk-end_trk).Unit();
407  const double tzen_trk = acos(dir_trk.Y()) * 180 / M_PI;
408  double tazi_trk = atan2(-dir_trk.X(), dir_trk.Z()) * 180 / M_PI + 332.066111;
409  if(tazi_trk < 0) tazi_trk += 360;
410  if(tazi_trk > 360) tazi_trk -= 360;
411 
412  const double tdotSun = AnglesToVector(zen, azi).Dot(AnglesToVector(tzen_trk, tazi_trk));
413  //const double tdotSunR = AnglesToVector(zenR, aziR).Dot(AnglesToVector(tzen_trk, tazi_trk));
414 
415  float track_entries[55] =
416  {
417  (float)e.id().run(), (float)e.id().subRun(), (float)e.id().event(),
418  slice, (float)i_track, (float)trackHits.size(), (float)sortedTrackHits.size(),
419  (float)outliers.size(), (float)P_up, (float)P_dn, (float)llr,
420  (float)chi2, (float)slope, (float)llrX, (float)chi2X, (float)slopeX,
421  (float)llrY, (float)chi2Y, (float)slopeY, (float)r2x, (float)r2y,
422  start.X(), start.Y(), start.Z(), start.T(), end.X(), end.Y(),
423  end.Z(), end.T(), (float)nxhit,(float)nyhit, dist, tDirX, tDirY,
424  tDirZ, tEleAngle, trackTotalE, containment, avgTX, avgTY,
425  (float)raw_daqheader->TotalMicroSlices(), (float)zen, (float)azi,
426  tCosTheta, tAzimuth, (float)tdotSun, (float)tzen_trk, (float)tazi_trk,
427  (float)event_time, (float)chi2X_fit, (float)chi2Y_fit,(float)trackX_fromfit, (float)trackY_fromfit, (float)trackX_fromfit_slice, (float)trackY_fromfit_slice
428  };
429  ntp_track->Fill(track_entries);
430  } //end loop track
431 }
float T() const
Definition: RecoHit.h:39
const art::PtrVector< rb::CellHit > & XCells() const
Get all cells from the x-view.
Definition: Cluster.h:124
float containmentType(rb::Track)
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
trk::CosmicTrackUtilities fUtilities
geo::View_t View() const
Definition: CellHit.h:41
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
float Z() const
Definition: RecoHit.h:38
A collection of associated CellHits.
Definition: Cluster.h:47
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
virtual TVector3 Start() const
Definition: Prong.h:73
T acos(T number)
Definition: d0nt_math.hpp:54
bool convertNovaTimeToUnixTime(uint64_t const &inputNovaTime, struct timespec &outputUnixTime)
double chi2()
#define M_PI
Definition: SbMath.h:34
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
Float_t Y
Definition: plot.C:38
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
double dist
Definition: runWimpSim.h:113
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
T atan(T number)
Definition: d0nt_math.hpp:66
Float_t Z
Definition: plot.C:38
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
unsigned short Cell() const
Definition: CellHit.h:40
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
double getLLR(std::set< std::pair< rb::CellHit, rb::RecoHit >, bool(*)(std::pair< rb::CellHit, rb::RecoHit >, std::pair< rb::CellHit, rb::RecoHit >)>, std::vector< rb::RecoHit > &outliers, double &P_up, double &P_dn, double &chi2, double &slope)
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
const art::PtrVector< rb::CellHit > & YCells() const
Get all cells from the x-view.
Definition: Cluster.h:126
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
float Azimuth(float const &dcosX, float const &dcosZ, float const &magnitude)
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
float CosTheta(float const &dcosY, float const &magnitude)
reference at(size_type n)
Definition: PtrVector.h:365
void GetSunPosition_FD(time_t time, double &zen, double &azi)
TVector3 Unit() const
const ana::Var wgt
size_type size() const
Definition: PtrVector.h:308
TH1F * h2
Definition: plot.C:45
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
TH1F * h1
bool pairHitComp(std::pair< rb::CellHit, rb::RecoHit > lhp, std::pair< rb::CellHit, rb::RecoHit > rhp)
float GeV() const
Definition: RecoHit.cxx:69
TVector3 AnglesToVector(double zen, double azi) const
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
Definition: MathUtil.cxx:36
virtual TVector3 InterpolateDir(double z) const
Definition: Track.cxx:348
double T
Definition: Xdiff_gwt.C:5
void LinFit(const std::vector< double > &x, const std::vector< double > &y, double *fitpar)
bool GoodTiming() const
Definition: CellHit.h:48
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Float_t e
Definition: plot.C:35
Float_t X
Definition: plot.C:38
art::ServiceHandle< locator::CelestialLocator > fSunPos
virtual bool Is3D() const
Definition: Prong.h:71
T atan2(T number)
Definition: d0nt_math.hpp:72
TVector3 upmuana::UpMuAnalysis::AnglesToVector ( double  zen,
double  azi 
) const
private

Definition at line 434 of file UpMuAnalysis_module.cc.

References std::cos(), M_PI, and std::sin().

Referenced by analyze().

435 {
436  zen *= M_PI / 180;
437  azi *= M_PI / 180;
438  // Doesn't matter which axis is which
439  return TVector3(sin(zen)*cos(azi), sin(zen)*sin(azi), cos(zen));
440 }
#define M_PI
Definition: SbMath.h:34
T sin(T number)
Definition: d0nt_math.hpp:132
T cos(T number)
Definition: d0nt_math.hpp:78
void upmuana::UpMuAnalysis::beginJob ( )
overridevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 144 of file UpMuAnalysis_module.cc.

References art::TFileDirectory::make(), and ntp_track.

145 {
146  // Implementation of optional member function here.
148  ntp_track = tfs->make<TNtuple>( "ntp_track", "Track Ntuple", "Run:SubRun:Event:SliceID:TrackID:Nhits:NRecohits:NOutliers:ProbUp:ProbDn:LLR:Chi2:Slope:LLRX:Chi2X:SlopeX:LLRY:Chi2Y:SlopeY:R2X:R2Y:StartX:StartY:StartZ:StartT:EndX:EndY:EndZ:EndT:TrackHitsX:TrackHitsY:Length:dirX:dirY:dirZ:eleAngle:totalE:containment:avgTX:avgTY:totalMSlices:SunZen:SunAzi:CosTheta:Azimuth:dotSun:zen_trk:azi_trk:event_time:chi2X_fit:chi2Y_fit:trackX_fromfit:trackY_fromfit:trackX_fromfit_slice:trackY_fromfit_slice");
149 
150 }
T * make(ARGS...args) const
detail::CachedProducts& art::EventObserverBase::cachedProducts ( )
inlineprotectedinherited

Definition at line 79 of file EventObserverBase.h.

References art::EventObserverBase::selectors_.

80  {
81  return selectors_;
82  }
detail::CachedProducts selectors_
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::consumes ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::consumes ( InputTag const &  it)
inherited

Definition at line 146 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

147 {
148  if (!moduleContext_)
149  return ProductToken<T>::invalid();
150 
151  consumables_[BT].emplace_back(ConsumableType::Product,
152  TypeID{typeid(T)},
153  it.label(),
154  it.instance(),
155  it.process());
156  return ProductToken<T>{it};
157 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 161 of file Consumer.h.

References T.

162 {
163  if (!moduleContext_)
164  return;
165 
166  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
167 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 171 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

172 {
173  if (!moduleContext_)
174  return ViewToken<T>::invalid();
175 
176  consumables_[BT].emplace_back(ConsumableType::ViewElement,
177  TypeID{typeid(T)},
178  it.label(),
179  it.instance(),
180  it.process());
181  return ViewToken<T>{it};
182 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
float upmuana::UpMuAnalysis::containmentType ( rb::Track  theTrack)
private

Definition at line 442 of file UpMuAnalysis_module.cc.

References fMaxX, fMaxY, fMaxZ, fMinX, fMinZ, febshutoff_auto::start, rb::Prong::Start(), and rb::Track::Stop().

Referenced by analyze().

443 {
444  TVector3 start = theTrack.Start();
445  TVector3 stop = theTrack.Stop();
446 
447  if (start.y() > stop.y()) // assume upward-going, so switch
448  { TVector3 hold = start; start = stop; stop = hold; }
449 
450  if (start.y() < fMinY || start.x() < fMinX || start.x() > fMaxX ||
451  start.z() < fMinZ || start.z() > fMaxZ) { // entered from outside
452  if (stop.y() > fMaxY ||
453  stop.x() < fMinX ||
454  stop.x() > fMaxX ||
455  stop.z() < fMinZ ||
456  stop.z() > fMaxZ) { return 1; }//through-going
457  return 2; // stopped
458  } // started in detector
459  if (stop.y() > fMaxY || stop.x() < fMinX || stop.x() > fMaxX ||
460  stop.z() < fMinZ || stop.z() > fMaxZ)
461  { return 3; }// in-produced
462 
463  return 4; // fully contained
464 
465 }
virtual TVector3 Start() const
Definition: Prong.h:73
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited
CurrentProcessingContext const* art::EDAnalyzer::currentContext ( ) const
protectedinherited
seed_t art::EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited
double upmuana::UpMuAnalysis::getLLR ( std::set< std::pair< rb::CellHit, rb::RecoHit >, bool(*)(std::pair< rb::CellHit, rb::RecoHit >, std::pair< rb::CellHit, rb::RecoHit >)>  ,
std::vector< rb::RecoHit > &  outliers,
double &  P_up,
double &  P_dn,
double &  chi2,
double &  slope 
)
private

Definition at line 467 of file UpMuAnalysis_module.cc.

References dist, upmuana::getDist(), upmuana::getErr(), in, LLR(), and test_ParserArtEvents::log.

Referenced by analyze().

474 {
475  if (sortedHits.size() < 2) {
476  P_up = 1e-30;
477  P_dn = 1e-30;
478  chi2 = 1e6;
479  slope = -1e6;
480  return 0;
481  }
482 
483  std::vector<double> measuredTimes;
484  std::vector<double> expectedTimes;
485  std::vector<double> wgts;
486 
487  // TGraph *g_plot = new TGraph (columns.size(), &columns[0], &columns[1]);
488  // g->SetMarkerStyle(21);
489  // g->Draw("alp");
490  // first hit
491  measuredTimes.push_back(0.0);
492  expectedTimes.push_back(0.0);
493  double err = getErr(sortedHits.begin()->first.PE());
494  wgts.push_back(1.0/err/err);
495 
496  double startX = sortedHits.begin()->second.X(),
497  startY = sortedHits.begin()->second.Y(),
498  startZ = sortedHits.begin()->second.Z(),
499  startT = sortedHits.begin()->second.T();
500 
501  std::vector<rb::RecoHit> in;
502 
503  for (auto i_hit = ++sortedHits.begin();
504  i_hit != sortedHits.end();
505  ++i_hit) {
506  in.push_back(i_hit->second);
507  measuredTimes.push_back(i_hit->second.T() - startT);
508  double dist = getDist(i_hit->second.X(), i_hit->second.Y(),
509  i_hit->second.Z(), startX, startY, startZ);
510  expectedTimes.push_back(dist/29.97);
511  err = getErr(i_hit->first.PE());
512  wgts.push_back(1.0/err/err);
513  }
514 
515  LLR(expectedTimes, measuredTimes, wgts, slope, chi2, P_up, P_dn,
516  in, outliers);
517 
518  return log(P_up/P_dn);
519 }
double chi2()
double getErr(double PE)
double dist
Definition: runWimpSim.h:113
void LLR(std::vector< double > &eT, std::vector< double > &mT, std::vector< double > &mTErr, double &slope, double &chi2, double &P_up, double &P_dn, std::vector< rb::RecoHit > &in, std::vector< rb::RecoHit > &outliers)
ifstream in
Definition: comparison.C:7
double getDist(double x1, double y1, double z1, double x2, double y2, double z2)
Float_t e
Definition: plot.C:35
art::Handle<art::TriggerResults> art::EventObserverBase::getTriggerResults ( Event const &  e) const
inlineinherited

Definition at line 61 of file EventObserverBase.h.

References art::detail::CachedProducts::getOneTriggerResults(), and art::EventObserverBase::selectors_.

62  {
64  }
detail::CachedProducts selectors_
art::Handle< art::TriggerResults > getOneTriggerResults(Event const &) const
Float_t e
Definition: plot.C:35
void upmuana::UpMuAnalysis::LinFit ( const std::vector< double > &  x,
const std::vector< double > &  y,
double *  fitpar 
)
private

Definition at line 615 of file UpMuAnalysis_module.cc.

References a, b, getGoodRuns4SAM::n, r(), std::sqrt(), submit_syst::x, xx, and submit_syst::y.

Referenced by analyze().

615  {
616 
617  const auto n = x_val.size();
618  const auto x = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/n; // <x>
619  const auto y = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/n; // <y>
620  const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/n; // <xx>
621  const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/n; // <yy>
622  const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/n; // <xy>
623 
624  const auto b = (xy - x*y)/(xx - x*x); // slope
625  const auto a = y - b*x; // intercept
626  const auto r = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y)); // Rxy - coeffcient of determination
627  fitpar[0] = a;
628  fitpar[1] = b;
629  fitpar[2] = r*r;
630 
631 }
Double_t xx
Definition: macro.C:12
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double a
const hit & b
Definition: hits.cxx:21
TRandom3 r(0)
void upmuana::UpMuAnalysis::LinFit ( const std::vector< double > &  x,
const std::vector< double > &  y,
const std::vector< double > &  ye,
double *  fitpar 
)
private

Definition at line 633 of file UpMuAnalysis_module.cc.

References a, b, DEFINE_ART_MODULE(), MECModelEnuComparisons::i, getGoodRuns4SAM::n, r(), std::sqrt(), submit_syst::x, xx, and submit_syst::y.

633  {
634 
635  int n = x_val.size();
636  double x = 0;
637  double y = 0;
638  double xx = 0;
639  double yy = 0;
640  double xy = 0;
641  double ee = 0;
642 
643  for ( int i=0; i<n; i++ ){
644  x = x + x_val[i]/y_err[i]/y_err[i];
645  y = y + y_val[i]/y_err[i]/y_err[i];
646 
647  xx = xx + x_val[i]*x_val[i]/y_err[i]/y_err[i];
648  yy = yy + y_val[i]*y_val[i]/y_err[i]/y_err[i];
649  xy = xy + x_val[i]*y_val[i]/y_err[i]/y_err[i];
650  ee = ee + 1./y_err[i]/y_err[i];
651  }
652 
653  const auto b = (xy*ee - x*y)/(xx*ee - x*x); // slope
654  const auto a = (xy - b*xx)/x; // intercept
655  const auto r = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y)); // Rxy - coeffcient of determination
656  fitpar[0] = a;
657  fitpar[1] = b;
658  fitpar[2] = r*r;
659 
660 
661 }
Double_t xx
Definition: macro.C:12
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double a
const hit & b
Definition: hits.cxx:21
TRandom3 r(0)
void upmuana::UpMuAnalysis::LLR ( std::vector< double > &  eT,
std::vector< double > &  mT,
std::vector< double > &  mTErr,
double &  slope,
double &  chi2,
double &  P_up,
double &  P_dn,
std::vector< rb::RecoHit > &  in,
std::vector< rb::RecoHit > &  outliers 
)
private

Definition at line 521 of file UpMuAnalysis_module.cc.

References a, b, stan::math::fabs(), stan::math::gamma_q(), MECModelEnuComparisons::i, util::LinFit(), getGoodRuns4SAM::n, cet::pow(), art::DataViewImpl::size(), and std::sqrt().

Referenced by getLLR().

527 {
528  // eT - x, mT - y
529  size_t n = mT.size();
530 
531  // y = a + bx
532  double a, b;
533  if (eT.size() >= 2)
534  chi2 = util::LinFit(eT, mT, wgts, b, a);
535  else {
536  chi2 = 1e6;
537  slope = -1e6;
538  P_up = 1e-30;
539  P_dn = 1e-30;
540  return;
541  }
542 
543  size_t totAllowOutlier = n/10;
544  size_t nOutliers = 0;
545  std::vector<double> x_filt, y_filt, w_filt;
546  for (size_t i=0; i < n; i++) {
547  double y_fit = a + b*eT.at(i);
548  if ( fabs(mT.at(i) - y_fit) < 5*(1.0/sqrt(wgts.at(i)))
549  || nOutliers>totAllowOutlier)
550  {
551  x_filt.push_back(eT.at(i));
552  y_filt.push_back(mT.at(i));
553  w_filt.push_back(wgts.at(i));
554  }
555  else {
556  ++nOutliers;
557  outliers.push_back(in[i]);
558  }
559  }
560 
561  if (x_filt.size() >= 2)
562  chi2 = util::LinFit(x_filt, y_filt, w_filt, b, a);
563  else {
564  chi2 = 1e6;
565  slope = -1e6;
566  P_up = 1e-30;
567  P_dn = 1e-30;
568  }
569  n = x_filt.size();
570  if (n < 5) {
571  slope = 0;
572  chi2 = 999;
573  P_up = 1e-30;
574  P_dn = 1;
575  return;
576  }
577 
578  slope = b;
579  chi2 /= (double)(n-2);
580 
581  double one_over_ee = 0.0,
582  x_over_ee = 0.0,
583  y_over_ee = 0.0;
584 
585  for (size_t i=0; i<n; ++i) {
586  one_over_ee += w_filt.at(i);
587  x_over_ee += x_filt.at(i)*w_filt.at(i);
588  y_over_ee += y_filt.at(i)*w_filt.at(i);
589  }
590 
591  double up_inter = (y_over_ee-x_over_ee)/one_over_ee;
592  double dn_inter = (y_over_ee+x_over_ee)/one_over_ee;
593 
594  double up_chi2 = 0.0, dn_chi2 = 0.0;
595  for (size_t i=0; i<n; ++i) {
596  double e = 1.0/sqrt(w_filt.at(i));
597  double up_expected = up_inter + x_filt.at(i);
598  double dn_expected = dn_inter - x_filt.at(i);
599  up_chi2 += pow((y_filt.at(i)-up_expected) / e, 2.0);
600  dn_chi2 += pow((y_filt.at(i)-dn_expected) / e, 2.0);
601  }
602 
603  double prob_up = boost::math::gamma_q((double)(n-2)/2.0,up_chi2/2.0),
604  prob_dn = boost::math::gamma_q((double)(n-2)/2.0,dn_chi2/2.0);
605 
606  if (prob_up < 1e-30) prob_up = 1e-30;
607  if (prob_dn < 1e-30) prob_dn = 1e-30;
608 
609  P_up = prob_up;
610  P_dn = prob_dn;
611 
612 }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
double chi2()
const double a
ifstream in
Definition: comparison.C:7
const hit & b
Definition: hits.cxx:21
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
Definition: MathUtil.cxx:36
Float_t e
Definition: plot.C:35
fvar< T > gamma_q(const fvar< T > &x1, const fvar< T > &x2)
Definition: gamma_q.hpp:12
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::mayConsume ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::mayConsume ( InputTag const &  it)
inherited

Definition at line 189 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

190 {
191  if (!moduleContext_)
192  return ProductToken<T>::invalid();
193 
194  consumables_[BT].emplace_back(ConsumableType::Product,
195  TypeID{typeid(T)},
196  it.label(),
197  it.instance(),
198  it.process());
199  return ProductToken<T>{it};
200 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 204 of file Consumer.h.

References T.

205 {
206  if (!moduleContext_)
207  return;
208 
209  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
210 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it)
inherited

Definition at line 214 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

215 {
216  if (!moduleContext_)
217  return ViewToken<T>::invalid();
218 
219  consumables_[BT].emplace_back(ConsumableType::ViewElement,
220  TypeID{typeid(T)},
221  it.label(),
222  it.instance(),
223  it.process());
224  return ViewToken<T>{it};
225 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
bool art::EventObserverBase::modifiesEvent ( ) const
inlineinherited

Definition at line 25 of file EventObserverBase.h.

26  {
27  return false;
28  }
static cet::exempt_ptr<Consumer> art::Consumer::non_module_context ( )
staticinherited
UpMuAnalysis& upmuana::UpMuAnalysis::operator= ( UpMuAnalysis const &  )
delete
UpMuAnalysis& upmuana::UpMuAnalysis::operator= ( UpMuAnalysis &&  )
delete
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited
std::string const& art::EventObserverBase::processName ( ) const
inlineinherited
void art::EventObserverBase::registerProducts ( MasterProductRegistry ,
ProductDescriptions ,
ModuleDescription const &   
)
inlineinherited

Definition at line 33 of file EventObserverBase.h.

References string.

36  {}
fhicl::ParameterSetID art::EventObserverBase::selectorConfig ( ) const
inlineinherited

Definition at line 56 of file EventObserverBase.h.

References art::EventObserverBase::selector_config_id_.

57  {
58  return selector_config_id_;
59  }
fhicl::ParameterSetID selector_config_id_
void art::Consumer::showMissingConsumes ( ) const
protectedinherited

Referenced by art::RootOutput::endJob().

void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited
bool art::EventObserverBase::wantAllEvents ( ) const
inlineinherited

Definition at line 46 of file EventObserverBase.h.

References art::EventObserverBase::wantAllEvents_.

Referenced by art::RootOutput::RootOutput().

47  {
48  return wantAllEvents_;
49  }
bool art::EventObserverBase::wantEvent ( Event const &  e)
inlineinherited

Definition at line 51 of file EventObserverBase.h.

References art::EventObserverBase::selectors_, and art::detail::CachedProducts::wantEvent().

52  {
53  return selectors_.wantEvent(e);
54  }
detail::CachedProducts selectors_
Float_t e
Definition: plot.C:35
bool wantEvent(Event const &)
std::string art::EDAnalyzer::workerType ( ) const
inlineinherited

Definition at line 109 of file EDAnalyzer.h.

References art::EDAnalyzer::currentContext().

110  {
111  return "WorkerT<EDAnalyzer>";
112  }

Member Data Documentation

double upmuana::UpMuAnalysis::fMaxX = 650.0
private

Definition at line 104 of file UpMuAnalysis_module.cc.

Referenced by containmentType().

double upmuana::UpMuAnalysis::fMaxY = 400.0
private

Definition at line 104 of file UpMuAnalysis_module.cc.

Referenced by containmentType().

double upmuana::UpMuAnalysis::fMaxZ = 5800.0
private

Definition at line 104 of file UpMuAnalysis_module.cc.

Referenced by containmentType().

double upmuana::UpMuAnalysis::fMinX = -650.0
private

Definition at line 103 of file UpMuAnalysis_module.cc.

Referenced by containmentType().

double upmuana::UpMuAnalysis::fMinY = -700.0
private

Definition at line 103 of file UpMuAnalysis_module.cc.

double upmuana::UpMuAnalysis::fMinZ = 100.0
private

Definition at line 103 of file UpMuAnalysis_module.cc.

Referenced by containmentType().

art::InputTag upmuana::UpMuAnalysis::fSliceInstance
private

Definition at line 100 of file UpMuAnalysis_module.cc.

art::InputTag upmuana::UpMuAnalysis::fSliceModuleLabel
private

Definition at line 99 of file UpMuAnalysis_module.cc.

art::ServiceHandle<locator::CelestialLocator> upmuana::UpMuAnalysis::fSunPos
private

Definition at line 126 of file UpMuAnalysis_module.cc.

Referenced by analyze().

trk::CosmicTrackUtilities upmuana::UpMuAnalysis::fUtilities
private

Definition at line 127 of file UpMuAnalysis_module.cc.

Referenced by analyze().

TNtuple* upmuana::UpMuAnalysis::ntp_track
private

Definition at line 129 of file UpMuAnalysis_module.cc.

Referenced by analyze(), and beginJob().

art::InputTag upmuana::UpMuAnalysis::trackTag_
private

Definition at line 98 of file UpMuAnalysis_module.cc.


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