20 #include "Utilities/func/MathUtil.h" 31 #include <boost/math/special_functions.hpp> 35 #include "NovaTimingUtilities/TimingUtilities.h" 46 #include "TGraphErrors.h" 51 return 165143/(1882.9+
pow(PE,2.11447)) + 10.4321;
55 double x2,
double y2,
double z2) {
56 return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
61 if (lhs.
Y() == rhs.
Y()) {
62 if (lhs.
Z() == rhs.
Z()) {
63 return lhs.
T() < rhs.
T();
65 else return lhs.
Z() < rhs.
Z();
67 else return lhs.
Y() < rhs.
Y();
71 std::pair<rb::CellHit, rb::RecoHit> rhp) {
106 double getLLR(std::set< std::pair<rb::CellHit,rb::RecoHit>,
107 bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
108 std::pair<rb::CellHit,rb::RecoHit>)> ,
109 std::vector<rb::RecoHit> &outliers,
110 double &P_up,
double &P_dn,
double &
chi2,
double &slope);
112 void LLR(std::vector<double>& eT,
113 std::vector<double>& mT,
114 std::vector<double>& mTErr,
double& slope,
double&
chi2,
115 double& P_up,
double& P_dn, std::vector<rb::RecoHit> &
in,
116 std::vector<rb::RecoHit> &outliers);
118 void LinFit(
const std::vector<double>&
x,
119 const std::vector<double>&
y,
double *fitpar);
120 void LinFit(
const std::vector<double>&
x,
121 const std::vector<double>&
y,
122 const std::vector<double>& ye,
double *fitpar);
139 , fSliceModuleLabel (p.
get<
std::
string>(
"SliceModuleLabel"))
140 , fSliceInstance (p.
get<
std::
string>(
"SliceInstanceLabel"))
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");
159 for(
int i = 0;
i < (
int)slicecol->size();
i++){
171 if (raw_triggers->size() != 1)
173 << __FILE__ <<
":" << __LINE__ <<
"\n" 174 <<
"RawTrigger vector has incorrect size: " << raw_triggers->size()
177 std::vector<const rb::Cluster*> slices;
178 std::vector<float> slice_containment;
184 unsigned long long event_time(raw_triggers->front().fTriggerTimingMarker_TimeStart);
195 std::vector<unsigned> tracks_per_slice;
196 for (
size_t i_track=0; i_track < tracks->size(); ++i_track){
198 rb::Track theTrack = tracks->at(i_track);
202 bool found_slice =
false;
205 const rb::Cluster *theSlice = &(*fmSlice.at(i_track));
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);
212 if(containment != 4) slice_containment.at(i_slice) = 1;
218 slice = (
float)slices.size();
219 slices.push_back(theSlice);
220 tracks_per_slice.push_back(1);
222 slice_containment.push_back(1);
224 slice_containment.push_back(4);
229 for(
size_t i_slice = 0; i_slice < slces.
size(); i_slice++)
231 if(slces.
at(i_slice).get() == theSlice)
233 slice = (
float)i_slice;
239 if (!theTrack.
Is3D())
continue;
242 std::set< std::pair<rb::CellHit,rb::RecoHit>,
243 bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
244 std::pair<rb::CellHit,rb::RecoHit>)>
246 std::set< std::pair<rb::CellHit,rb::RecoHit>,
247 bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
248 std::pair<rb::CellHit,rb::RecoHit>)>
250 std::set< std::pair<rb::CellHit,rb::RecoHit>,
251 bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
252 std::pair<rb::CellHit,rb::RecoHit>)>
255 std::vector<double> x_hit;
256 std::vector<double> y_hit;
257 std::vector<double> zy_hit;
258 std::vector<double> zx_hit;
261 for (
size_t i_hit=0; i_hit < trackHits.
size(); ++i_hit) {
268 unsigned short iC = theHit->
Cell();
269 unsigned short iP = theHit->
Plane();
273 zx_hit.push_back(iP);
278 zy_hit.push_back(iP);
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();
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;
300 Double_t trackX_fromfit_slice = 0;
301 Double_t trackY_fromfit_slice = 0;
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);
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;
313 std::vector<double> weight_x;
314 std::vector<double> weight_y;
322 for(
size_t y_ind =0; y_ind < ycells.
size(); y_ind++)
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));
335 for(
size_t x_ind = 0; x_ind < xcells.
size(); x_ind++)
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));
347 trackY_fromfit_slice += trackY_fromfit;
348 trackX_fromfit_slice += trackX_fromfit;
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;
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;
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;
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);
371 if (sortedTrackHits.size() >= 1) start = (sortedTrackHits.begin()->second);
373 if (sortedTrackHits.size() >= 2) end = (--sortedTrackHits.end())->
second;
378 float tDirX=0, tDirY=0, tDirZ=0, tEleAngle, trackTotalE=0;
379 float avgTX=0, avgTY=0, tCosTheta, tAzimuth;
381 for (
size_t i_hit=0; i_hit < trackHits.
size(); ++i_hit) {
387 tDirX += theDir.x()/trackHits.
size();
388 tDirY += theDir.y()/trackHits.
size();
389 tDirZ += theDir.z()/trackHits.
size();
394 avgTX += theRecoHit.
T()/sortedTrackHitsX.size();
396 avgTY += theRecoHit.
T()/sortedTrackHitsY.size();
399 tEleAngle =
atan(tDirY/
sqrt(tDirX*tDirX + tDirZ*tDirZ));
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;
415 float track_entries[55] =
418 slice, (
float)i_track, (
float)trackHits.
size(), (
float)sortedTrackHits.size(),
420 (
float)
chi2, (
float)slope, (
float)llrX, (
float)chi2X, (
float)slopeX,
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,
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
445 TVector3 stop = theTrack.
Stop();
447 if (start.y() > stop.y())
448 { TVector3 hold =
start; start = stop; stop = hold; }
450 if (start.y() < fMinY || start.x() <
fMinX || start.x() >
fMaxX ||
452 if (stop.y() >
fMaxY ||
456 stop.z() >
fMaxZ) {
return 1; }
468 bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
469 std::pair<rb::CellHit,rb::RecoHit>)>
471 std::vector<rb::RecoHit> &outliers,
472 double &P_up,
double &P_dn,
double &
chi2,
475 if (sortedHits.size() < 2) {
483 std::vector<double> measuredTimes;
484 std::vector<double> expectedTimes;
485 std::vector<double> wgts;
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);
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();
501 std::vector<rb::RecoHit>
in;
503 for (
auto i_hit = ++sortedHits.begin();
504 i_hit != sortedHits.end();
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);
515 LLR(expectedTimes, measuredTimes, wgts, slope, chi2, P_up, P_dn,
518 return log(P_up/P_dn);
522 std::vector<double>& mT,
523 std::vector<double>& wgts,
double& slope,
524 double&
chi2,
double& P_up,
double& P_dn,
525 std::vector<rb::RecoHit>&
in,
526 std::vector<rb::RecoHit>& outliers)
529 size_t n = mT.size();
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)
551 x_filt.push_back(eT.at(
i));
552 y_filt.push_back(mT.at(
i));
553 w_filt.push_back(wgts.at(
i));
557 outliers.push_back(in[
i]);
561 if (x_filt.size() >= 2)
579 chi2 /= (double)(n-2);
581 double one_over_ee = 0.0,
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);
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;
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);
606 if (prob_up < 1e-30) prob_up = 1e-30;
607 if (prob_dn < 1e-30) prob_dn = 1e-30;
617 const auto n = x_val.size();
618 const auto x = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/
n;
619 const auto y = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/
n;
620 const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/
n;
621 const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/
n;
622 const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/
n;
624 const auto b = (xy -
x*
y)/(
xx -
x*
x);
625 const auto a =
y -
b*
x;
626 const auto r = (xy - x*
y)/
sqrt((
xx - x*x)*(yy -
y*
y));
635 int n = x_val.size();
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];
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];
653 const auto b = (xy*ee - x*
y)/(xx*ee - x*x);
654 const auto a = (xy -
b*
xx)/x;
655 const auto r = (xy - x*
y)/
sqrt((xx - x*x)*(yy - y*
y));
const art::PtrVector< rb::CellHit > & XCells() const
Get all cells from the x-view.
float containmentType(rb::Track)
fvar< T > fabs(const fvar< T > &x)
Float_t y1[n_points_granero]
UpMuAnalysis & operator=(UpMuAnalysis const &)=delete
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
trk::CosmicTrackUtilities fUtilities
Float_t x1[n_points_granero]
Vertical planes which measure X.
art::InputTag fSliceModuleLabel
A collection of associated CellHits.
A rb::Prong with full reconstructed trajectory.
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
bool convertNovaTimeToUnixTime(uint64_t const &inputNovaTime, struct timespec &outputUnixTime)
Horizontal planes which measure Y.
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
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 swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
Locate positions of celestial bodies.
virtual double TotalLength() const
Length (cm) of all the track segments.
UpMuAnalysis(fhicl::ParameterSet const &p)
unsigned short Cell() const
void push_back(Ptr< U > const &p)
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)
bool recoHitComp(rb::RecoHit lhs, rb::RecoHit rhs)
const art::PtrVector< rb::CellHit > & YCells() const
Get all cells from the x-view.
virtual TVector3 Dir() const
Unit vector describing prong direction.
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
float Azimuth(float const &dcosX, float const &dcosZ, float const &magnitude)
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
float CosTheta(float const &dcosY, float const &magnitude)
reference at(size_type n)
EDAnalyzer(Table< Config > const &config)
void GetSunPosition_FD(time_t time, double &zen, double &azi)
Vertex location in position and time.
EventNumber_t event() const
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
bool pairHitComp(std::pair< rb::CellHit, rb::RecoHit > lhp, std::pair< rb::CellHit, rb::RecoHit > rhp)
SubRunNumber_t subRun() const
T * make(ARGS...args) const
TVector3 AnglesToVector(double zen, double azi) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
void analyze(art::Event const &e) override
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...
virtual TVector3 InterpolateDir(double z) const
void LinFit(const std::vector< double > &x, const std::vector< double > &y, double *fitpar)
TVector3 Stop() const
Position of the final trajectory point.
double getDist(double x1, double y1, double z1, double x2, double y2, double z2)
art::ServiceHandle< locator::CelestialLocator > fSunPos
fvar< T > gamma_q(const fvar< T > &x1, const fvar< T > &x2)
virtual bool Is3D() const
art::InputTag fSliceInstance