21 #include "TTimeStamp.h" 37 #include <TGraphErrors.h> 38 #include <TPaveText.h> 42 #include <TApplication.h> 51 #include <boost/math/special_functions/gamma.hpp> 52 #include <boost/math/special_functions.hpp> 104 double GetT(
double tdc,
double dist_to_readout);
109 void LinFit(
const std::vector<double>&
x,
const std::vector<double>&
y,
double *fitpar);
110 void LinFit(
const std::vector<double>&
x,
const std::vector<double>&
y,
const std::vector<double>& ye,
double *fitpar);
111 void LinFitLLR(std::vector<double>& x_hit, std::vector<double>& y_hit, std::vector<double>& y_error,
112 double& slope,
double&
chi2,
double& P_up,
double& P_dn);
162 produces<std::vector<TriggerDecision>>();
172 std::unique_ptr<std::vector<TriggerDecision>>
173 trigger_decisions(
new std::vector<TriggerDecision>);
182 assert(fohl.size() == tracks->size());
184 for(
size_t i_track = 0; i_track < tracks->size(); ++i_track){
191 uint64_t min_tdc = this_hit_list->front().TDC().val;
192 uint64_t max_tdc = this_hit_list->back().TDC().val;
193 for(
size_t i_hit = 0; i_hit < this_hit_list->size(); ++i_hit){
194 if (this_hit_list->at(i_hit).TDC().val < min_tdc) min_tdc = this_hit_list->at(i_hit).TDC().val;
195 if (this_hit_list->at(i_hit).TDC().val > max_tdc) max_tdc = this_hit_list->at(i_hit).TDC().val;
200 int StartX = track.
Start().X();
201 int StartY = track.
Start().Y();
202 int StartZ = track.
Start().Z();
204 int EndX = track.
End().X();
205 int EndY = track.
End().Y();
206 int EndZ = track.
End().Z();
214 length.SetXYZ(
cw*(track.
End().X() - track.
Start().X()),
219 double tr_length = length.Mag();
220 double exp_trav_time = tr_length/29.97;
223 if (TMath::Abs(EndX - StartX) <
_dX)
continue;
224 if (TMath::Abs(EndY - StartY) <
_dY)
continue;
225 if (TMath::Abs(EndZ - StartZ) <
_dZ)
continue;
226 if (!track.
Is3D())
continue;
229 std::vector<double> x_hit;
230 std::vector<double> y_hit;
231 std::vector<double> zy_hit;
232 std::vector<double> zx_hit;
235 for(
size_t i_hit = 0; i_hit < this_hit_list->size(); ++i_hit){
236 uint16_t iC = this_hit_list->at(i_hit).Cell().val;
237 uint16_t iP = this_hit_list->at(i_hit).Plane().val;
238 uint8_t
view = this_hit_list->at(i_hit).View().val;
239 int iTDC =
static_cast<int>(this_hit_list->at(i_hit).TDC().val - absolute_t0.
val);
243 zx_hit.push_back(iP);
250 zy_hit.push_back(iP);
259 LinFit(x_hit, zx_hit, fitpar);
260 double r2x = fitpar[2];
261 LinFit(y_hit, zy_hit, fitpar);
262 double r2y = fitpar[2];
264 if(r2x <
_R2X)
continue;
265 if(r2y <
_R2Y)
continue;
268 std::vector<double> eT_vecxy;
269 std::vector<double> mT_vecxy;
270 std::vector<double> mT_vecx;
271 std::vector<double> mT_vecy;
274 std::vector<double> mT_vecxy_e;
276 for(
size_t i_hit = 0; i_hit < this_hit_list->size(); ++i_hit){
277 uint16_t iC = this_hit_list->at(i_hit).Cell().val;
278 uint16_t iP = this_hit_list->at(i_hit).Plane().val;
279 int hit_adc = this_hit_list->at(i_hit).ADC().val;
281 if (hit_adc < 50)
continue;
284 DAQHit const& dqhit = this_hit_list->at(i_hit);
287 TVector3
end = track.
End();
288 GetXYZ(dqhit, start, end, pos);
290 int iTDC =
static_cast<int>(this_hit_list->at(i_hit).TDC().val - absolute_t0.
val);
291 double meas_time =
static_cast<double>(iTDC -
T0);
292 meas_time = meas_time + 0.01*
static_cast<double>((this_hit_list->at(i_hit).TDC().frac));
294 double tns =
GetT(meas_time, pos[3] + pigtail);
296 double DCM_off =
getDCMOff(iC, this_hit_list->at(i_hit).View().val);
300 if(iP > TMath::Max(StartZ, EndZ))
continue;
301 if(iP < TMath::Min(StartZ, EndZ))
continue;
303 if(iC > TMath::Max(StartX, EndX))
continue;
304 if(iC < TMath::Min(StartX, EndX))
continue;
307 if(iP > TMath::Max(StartZ, EndZ))
continue;
308 if(iP < TMath::Min(StartZ, EndZ))
continue;
310 if(iC > TMath::Max(StartY, EndY))
continue;
311 if(iC < TMath::Min(StartY, EndY))
continue;
314 double exp_time = -999;
319 exp_time = exp_trav_time*(iC-StartX)/(EndX-StartX);
322 exp_time = exp_trav_time*(iC-StartY)/(EndY-StartY);
324 if(
false && tr_length > 1200. && tns<-1000){
326 std::cout <<
", i_track: " << std::setw(3) << i_track
327 <<
", i_hit: " << std::setw(3) << i_hit
328 <<
", iP: " << std::setw(3) << iP
329 <<
", iC: " << std::setw(3) << iC
330 <<
", T0: " << std::setw(3) << T0
331 <<
", pigtail: " << std::setw(3) << pigtail
332 <<
", DCM_off: " << std::setw(3) << DCM_off
333 <<
", ADC: " << std::setw(5) << (
int)this_hit_list->at(i_hit).ADC().val
334 <<
", TDC: " << (uint64_t)this_hit_list->at(i_hit).TDC().val
335 <<
", meas_time: " << meas_time
337 <<
", exp_time: " << exp_time
343 eT_vecxy.push_back(exp_time);
344 mT_vecxy.push_back(tns);
345 mT_vecxy_e.push_back(
GetRes(hit_adc));
347 mT_vecx.push_back(tns);
349 mT_vecy.push_back(tns);
353 unsigned nxhit = mT_vecx.size();
354 unsigned nyhit = mT_vecy.size();
361 double slope_xy, chi_xy, P_up_xy, P_dn_xy;
362 LinFitLLR(eT_vecxy, mT_vecxy, mT_vecxy_e, slope_xy, chi_xy, P_up_xy, P_dn_xy);
363 double LLR =
log(P_up_xy/P_dn_xy);
367 if(LLR <
_LLR)
continue;
369 if(chi_xy >
_Chi2)
continue;
376 std::cout<< std::setiosflags(std::ios::fixed)
377 <<
"Run: " << std::setprecision(4) << std::setw(8) <<
std::left << e.
id().
run()
378 <<
", Event: " << std::setprecision(4) << std::setw(8) <<
std::left << e.
id().
event()
379 <<
", LLR: " << std::setprecision(4) << std::setw(5) <<
std::left << LLR
380 <<
", Chi2: " << std::setprecision(4) << std::setw(5) <<
std::left << chi_xy
381 <<
", slope: " << std::setprecision(4) << std::setw(5) <<
std::left << slope_xy
382 <<
", tr_length: " << std::setprecision(3) << std::setw(3) <<
std::left << tr_length
383 <<
", dX: " << std::setprecision(4) << std::setw(4) <<
std::left << TMath::Abs(EndX - StartX)
384 <<
", dY: " << std::setprecision(4) << std::setw(4) <<
std::left << TMath::Abs(EndY - StartY)
385 <<
", dZ: " << std::setprecision(4) << std::setw(4) <<
std::left << TMath::Abs(EndZ - StartZ)
386 <<
", 3d: " << std::setprecision(1) << std::setw(1) <<
std::left << track.
Is3D()
387 <<
", nxhit: " << std::setprecision(4) << std::setw(3) <<
std::left << nxhit
388 <<
", nyhit: " << std::setprecision(4) << std::setw(3) <<
std::left << nyhit
389 <<
", nxyhit: " << std::setprecision(4) << std::setw(3) <<
std::left << nxhit+nyhit
390 <<
", r2x: " << std::setprecision(4) << std::setw(5) <<
std::left << r2x
391 <<
", r2y: " << std::setprecision(4) << std::setw(5) <<
std::left << r2y
409 trigger_decisions->emplace_back
412 trigger_decisions->emplace_back
421 e.
put(std::move(trigger_decisions));
449 double x0 =
c0x +
cw*start.X();
double x1 =
c0x +
cw*end.X();
450 double y0 =
c0y +
cw*start.Y();
double y1 =
c0y +
cw*end.Y();
451 double z0 =
p0 +
pw*start.Z();
double z1 =
p0 +
pw*end.Z();
454 y = (y1 - y0)/(z1 - z0);
460 x = (x1 - x0)/(z1 - z0);
475 double speedOfFiberTransport = 15.3;
476 double TDC_to_ns = 15.625;
481 return (tdc*TDC_to_ns - dist_to_readout/speedOfFiberTransport);
494 double res = p0/(p1 +
pow( static_cast<double>(ADC), p2) ) + p3;
513 double res = p0/(p1 +
pow( static_cast<double>(ADC), p2) ) + p3;
519 UInt_t cellid = daqhit.
Cell().
val;
523 const int kCellsPerModule = 32;
524 cellid = cellid % kCellsPerModule;
528 cellid = kCellsPerModule-cellid-1;
531 if(cellid >= kCellsPerModule)
return 100;
535 const double kPigtails[kCellsPerModule] = {
536 34.5738, 38.4379, 42.3020, 46.1660, 50.0301, 53.8942, 57.7582, 61.6223,
537 64.7504, 68.6144, 72.4785, 76.3426, 80.2067, 84.0707, 87.9348, 91.0790,
538 95.3301, 99.1941, 103.058, 106.922, 110.786, 114.650, 118.514, 122.379,
539 125.507, 129.371, 133.235, 137.099, 140.963, 144.827, 148.691, 150.751
541 return kPigtails[cellid];
550 const auto n = x_val.size();
551 const auto x = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/
n;
552 const auto y = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/
n;
553 const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/
n;
554 const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/
n;
555 const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/
n;
557 const auto b = (xy -
x*
y)/(
xx -
x*
x);
558 const auto a =
y -
b*
x;
559 const auto r = (xy - x*
y)/
sqrt((
xx - x*x)*(yy -
y*
y));
566 void novaddt::UpMuTrigger::LinFit(
const std::vector<double>& x_val,
const std::vector<double>& y_val,
const std::vector<double>& y_err,
double *fitpar){
568 int n = x_val.size();
576 for (
int i=0;
i<
n;
i++ ){
578 x = x + x_val[
i]/y_err[
i]/y_err[
i];
579 y = y + y_val[
i]/y_err[
i]/y_err[
i];
581 xx = xx + x_val[
i]*x_val[
i]/y_err[
i]/y_err[
i];
582 yy = yy + y_val[
i]*y_val[
i]/y_err[
i]/y_err[
i];
583 xy = xy + x_val[
i]*y_val[
i]/y_err[
i]/y_err[
i];
584 ee = ee + 1./y_err[
i]/y_err[
i];
587 const auto b = (xy*ee - x*
y)/(xx*ee - x*x);
588 const auto a = (xy -
b*
xx)/x;
589 const auto r = (xy - x*
y)/
sqrt((xx - x*x)*(yy - y*
y));
598 double& slope,
double&
chi2,
double& P_up,
double& P_dn) {
600 int n = x_hit.size();
603 LinFit(x_hit, y_hit, y_error, fitpar);
604 double a = fitpar[0];
605 double b = fitpar[1];
608 int totAllowOutlier =
static_cast<int>(0.1*
n);
609 const double sigma_cut = 5;
611 std::vector<double> x_filt, y_filt, ye_filt, x_filt_e;
612 for (
int i = 0;
i <
n;
i++) {
613 double y_fit = a + b*x_hit[
i];
614 if (
fabs(y_hit[
i] - y_fit) < sigma_cut*y_error[
i] || iOutlier>totAllowOutlier)
616 x_filt.push_back(x_hit[
i]);
617 y_filt.push_back(y_hit[i]);
618 ye_filt.push_back(y_error[i]);
619 x_filt_e.push_back(0);
626 LinFit(x_filt, y_filt, ye_filt, fitpar);
643 for (
int i=0;
i<
n;
i++) {
644 double y_expected = a + b*x_filt.at(
i);
645 chi2+=
pow((y_filt.at(
i)-y_expected) / ye_filt.at(
i), 2.0);
648 chi2/=
static_cast<double>(n-2);
652 double one_over_ee = 0.0;
653 double x_over_ee = 0.0;
654 double y_over_ee = 0.0;
656 for (
int i=0;
i<
n;
i++) {
657 double e = ye_filt.at(
i);
658 one_over_ee+=1.0/e/
e;
659 x_over_ee+=x_filt.at(
i)/e/
e;
660 y_over_ee+=y_filt.at(
i)/e/
e;
664 double up_inter = (y_over_ee-x_over_ee)/one_over_ee;
665 double down_inter = (y_over_ee+x_over_ee)/one_over_ee;
669 double up_chi2=0.0, down_chi2=0.0;
670 for (
int i=0;
i<
n;
i++) {
672 double e = ye_filt.at(
i);
673 double up_expected = up_inter + x_filt.at(
i);
674 double down_expected = down_inter - x_filt.at(
i);
675 up_chi2 +=
pow((y_filt.at(
i)-up_expected) / e, 2.0);
676 down_chi2 +=
pow((y_filt.at(
i)-down_expected) / e, 2.0);
688 if (prob_up<1
e-30) prob_up = 1
e-30;
689 if (prob_down<1
e-30) prob_down = 1
e-30;
696 TCanvas *
can =
new TCanvas(
"can",
"can", 800, 600);
698 TGraphErrors *eTmT_grx =
new TGraphErrors(x_filt.size(), &x_filt[0], &y_filt[0], &x_filt_e[0], &ye_filt[0]);
699 TF1 *SlopeFit =
new TF1(
"SlopeFit",
"pol1", -10, 150);
700 TF1 *UpFit =
new TF1(
"UpFit",
"pol1", -10, 150);
701 TF1 *UpFit1 =
new TF1(
"UpFit1",
"pol1", -10, 150);
702 TF1 *DnFit2 =
new TF1(
"DnFit2",
"pol1", -10, 150);
703 TF1 *DnFit =
new TF1(
"DnFit",
"pol1", -10, 150);
704 UpFit->SetLineColor(4);
705 DnFit->SetLineColor(4);
706 DnFit->SetLineStyle(3);
707 UpFit->SetParameter(0, up_inter);
708 UpFit->SetParameter(1, 1);
710 DnFit2->SetParameter(0, fitpar[0]);
711 DnFit2->SetParameter(1, fitpar[1]);
712 DnFit2->SetLineColor(
kGreen);
713 DnFit2->SetLineStyle(3);
716 DnFit->SetParameter(0, down_inter);
717 DnFit->SetParameter(1, -1);
721 SlopeFit->SetParameter(0, a);
722 SlopeFit->SetParameter(1, b);
723 SlopeFit->SetLineColor(3);
724 SlopeFit->SetLineStyle(3);
728 eTmT_grx->GetXaxis()->SetRangeUser(-10, 100);
729 eTmT_grx->Draw(
"AP");
730 SlopeFit->Draw(
"same");
732 DnFit2->Draw(
"same");
735 UpFit1->SetLineColor(
kRed);
736 UpFit1->SetLineWidth(1);
739 eTmT_grx->Fit(
"UpFit1",
"BR",
"", -10, 150);
741 TPaveText *
pt =
new TPaveText(.75, .9, 0.9, 1.0,
"brNDC");
742 pt->SetBorderSize(1);
743 pt->AddText(Form(
"P_up = %.1e", P_up ));
744 pt->AddText(Form(
"P_dn = %.1e", P_dn ));
745 pt->AddText(Form(
"LLR = %.1f",
log(P_up/P_dn) ));
746 pt->AddText(Form(
"Chi2 = %.2f", chi2 ));
747 pt->AddText(Form(
"NDF = %d", n-2 ));
752 can->SaveAs(Form(
"plots/Chi2_%f.pdf", chi2));
768 if (cell > 31 && cell < 64)
return -26;
769 else if (cell > 63)
return -40;
772 if (cell > 7 && cell < 32)
return -10;
773 else if (cell > 31)
return 5;
777 return (cell/64)*40.;
double GetT(double tdc, double dist_to_readout)
fvar< T > fabs(const fvar< T > &x)
Float_t y1[n_points_granero]
double getDCMOff(uint16_t cell, uint8_t view)
novaddt::Plane const & Plane() const
Float_t x1[n_points_granero]
void LinFit(const std::vector< double > &x, const std::vector< double > &y, double *fitpar)
DEFINE_ART_MODULE(TestTMapFile)
std::string _trackToSlice
TVector3 const & End() const
void LinFitLLR(std::vector< double > &x_hit, std::vector< double > &y_hit, std::vector< double > &y_error, double &slope, double &chi2, double &P_up, double &P_dn)
Identifier for the Y measuring view of the detector (side)
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
ProductID put(std::unique_ptr< PROD > &&product)
novaddt::View const & View() const
Identifier for the X measuring view of the detector (top)
TVector3 const & Start() const
std::string _trackinstance
bool const & Is3D() const
EventNumber_t event() const
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
UpMuTrigger(fhicl::ParameterSet const &p)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
virtual bool filter(art::Event &e) override
novaddt::Cell const & Cell() const
assert(nhit_max >=nhit_nbins)
std::string _sliceinstance
double GetPigtail(DAQHit const &daqhit) const
void GetXYZ(DAQHit const &daqhit, TVector3 start, TVector3 end, double *xyzd)
fvar< T > gamma_q(const fvar< T > &x1, const fvar< T > &x2)