21 #include "Utilities/func/MathUtil.h" 30 #include "NovaDAQConventions/DAQConventions.h" 49 class FmmTrackerValidation;
71 TRandom *
genR =
new TRandom3();
94 double doca2(
double x,
double y,
double x0,
double y0,
double dydx);
97 double vfit(std::vector<TVector3> & Positions,
98 std::vector<double> & Ts,
99 std::vector<double> & Terr,
104 void LinFit(
const std::vector<double>&
x,
const std::vector<double>&
y,
double *fitpar);
126 RCV_Tree = tfs->
make<TTree>(
"RCV_Tree",
"Information from RC track info");
253 if (p->
PdgCode()!=42 )
continue;
260 cosx = (p->
Px())/momentum;
261 cosy = (p->
Py())/momentum;
262 cosz = (p->
Pz())/momentum;
264 double momentum2 = (p->
P())*(p->
P());
265 double mass2 = (p->
Mass())*(p->
Mass());
268 unsigned ntruehits_Tot = flshits.size();
269 if ( ntruehits_Tot==0 )
break;
270 ti_MC = flshits[0].GetEntryT();
271 xi_MC = flshits[0].GetEntryX();
272 yi_MC = flshits[0].GetEntryY();
273 zi_MC = flshits[0].GetEntryZ();
274 tf_MC = flshits[0].GetExitT();
275 xf_MC = flshits[0].GetExitX();
276 yf_MC = flshits[0].GetExitY();
277 zf_MC = flshits[0].GetExitZ();
278 Edep_MC = flshits[0].GetEdep();
280 trlMC = flshits[0].GetTotalPathLength();
282 unsigned int planeID_i = flshits[0].GetPlaneID();
283 unsigned int cellID_i = flshits[0].GetCellID();
284 unsigned int planeID_f = planeID_i;
285 unsigned int cellID_f = cellID_i;
290 for (
unsigned h =1;
h!= flshits.size(); ++
h) {
291 double ti = flshits[
h].GetEntryT();
292 double tf = flshits[
h].GetExitT();
295 xi_MC = flshits[
h].GetEntryX();
296 yi_MC = flshits[
h].GetEntryY();
297 zi_MC = flshits[
h].GetEntryZ();
298 planeID_i = flshits[
h].GetPlaneID();
299 cellID_i = flshits[
h].GetCellID();
303 xf_MC = flshits[
h].GetExitX();
304 yf_MC = flshits[
h].GetExitY();
305 zf_MC = flshits[
h].GetExitZ();
306 planeID_f = flshits[
h].GetPlaneID();
307 cellID_f = flshits[
h].GetCellID();
314 trlMC += flshits[
h].GetTotalPathLength();
319 double local[3], world[3];
371 int max_nhits_true = 0;
372 int trakIdx_for_max_nhits_true = -1;
373 for (
unsigned trkIdx = 0; trkIdx != trks->size(); ++trkIdx ) {
377 for (
unsigned i = 0;
i != track.
NCell(); ++
i) {
383 trakIdx_for_max_nhits_true=trkIdx;
388 if (trakIdx_for_max_nhits_true>-1) {
392 for (
unsigned i = 0;
i != track.
NCell(); ++
i) {
395 double tns = chit->
TNS();
396 double adc = chit->
ADC();
400 if (track.
NCell()>=6){
406 for (
unsigned i = 0;
i != track.
NCell(); ++
i) {
408 double tns = chit->
TNS();
409 double adc = chit->
ADC();
423 std::vector<int> planes;
430 for (
unsigned i = 0;
i != track.
NCell(); ++
i) {
432 double tns = chit->
TNS();
436 planes.emplace_back(chit->
Plane());
441 if (chit->
Plane() < P_x0) P_x0 = chit->
Plane();
442 if (chit->
Plane() > P_x1) P_x1 = chit->
Plane();
446 if (chit->
Plane() < P_y0) P_y0 = chit->
Plane();
447 if (chit->
Plane() > P_y1) P_y1 = chit->
Plane();
451 if (planes.size()>=3 && P_x1-P_x0>=3 && P_y1-P_y0>=3) {
454 std::sort (planes.begin(), planes.end());
457 p_cross = planes.back()-planes.front()+1;
459 bool start_flag =
false;
460 for (
unsigned pidx =0; pidx+1<planes.size(); ++pidx) {
461 int dP = planes[pidx+1] - planes[pidx] - 1;
464 if ( (planes[pidx+1] - planes[pidx]) % 2 ) start_flag =
true;
467 if ( (planes[pidx+1] - planes[pidx]) % 2 ) ++
p_overlap;
471 if (dP > max_gap) max_gap = dP;
476 std::vector<double> x0s_x, x1s_x, y0s_y, y1s_y;
477 std::vector<double> x0s_z, x1s_z, y0s_z, y1s_z;
478 for (
unsigned i = 0;
i != track.
NCell(); ++
i) {
480 double tns = chit->
TNS();
485 if (chit->
Plane() == P_x0 ) {
486 x0s_x.emplace_back(xyz[0]);
487 x0s_z.emplace_back(xyz[2]);
489 if (chit->
Plane() == P_x1 ) {
490 x1s_x.emplace_back(xyz[0]);
491 x1s_z.emplace_back(xyz[2]);
494 if (chit->
Plane() == P_y0 ) {
495 y0s_y.emplace_back(xyz[1]);
496 y0s_z.emplace_back(xyz[2]);
498 if (chit->
Plane() == P_y1 ) {
499 y1s_y.emplace_back(xyz[1]);
500 y1s_z.emplace_back(xyz[2]);
514 for (
unsigned h0 = 0; h0 != x0s_x.size(); ++h0) {
515 for (
unsigned h1 = 0;
h1 != x1s_x.size(); ++
h1) {
517 dx = x1s_x[
h1]-x0s_x[h0];
518 dxdz = dx/(x1s_z[
h1]-x0s_z[h0]);
525 for (
unsigned h0 = 0; h0 != y0s_y.size(); ++h0) {
526 for (
unsigned h1 = 0;
h1 != y1s_y.size(); ++
h1) {
528 dy = y1s_y[
h1]-y0s_y[h0];
529 dydz = dy/(y1s_z[
h1]-y0s_z[h0]);
539 double xi = dxdz*(zi-z_x0)+x_x0;
540 double xf = dxdz*(zf-z_x0)+x_x0;
541 double yi = dydz*(zi-z_y0)+y_y0;
542 double yf = dydz*(zf-z_y0)+y_y0;
543 TVector3 Pi(xi,yi,zi);
544 TVector3 Pf(xf,yf,zf);
550 double t_min = 9999999;
551 double t_max = -9999999;
553 std::vector<double> ts_sat,terrs_sat;
554 std::vector<double> ts_unsat,terrs_unsat;
555 std::vector<TVector3> pos_sat, pos_unsat;
561 TVector3 adc_Center(0,0,0);
562 TVector3 GeV_Center(0,0,0);
563 TVector3 geo_Center(0,0,0);
566 for (
unsigned i = 0;
i != track.
NCell(); ++
i) {
568 double tns = chit->
TNS();
575 W = dydz*(xyz[2]-z_y0)+y_y0;
576 lxsqr +=
doca2(xyz[2],xyz[0],z_x0,x_x0, dxdz);
579 W = dxdz*(xyz[2]-z_x0)+x_x0;
580 lysqr +=
doca2(xyz[2],xyz[1],z_y0,y_y0, dydz);
584 if (!rhit.IsCalibrated())
continue;
585 float hitE = rhit.
GeV();
587 float hitT = rhit.T();
588 unsigned int hitCell = chit->
Cell();
589 unsigned int hitPlane = chit->
Plane();
591 float hitX = rhit.X();
594 hitxe.push_back(hitE);
595 hitxt.push_back(hitT);
596 hitxx.push_back(hitX);
599 float hitY = rhit.Y();
602 hitye.push_back(hitE);
603 hityt.push_back(hitT);
604 hityy.push_back(hitY);
610 TVector3
pos(rhit.X(), rhit.Y(), rhit.Z());
613 adc_Center += (chit->
ADC())*
pos;
614 GeV_Center += rhit.GeV() *
pos;
620 if (chit->
ADC(3) == 4095) {
621 pos_sat.emplace_back(
pos-Pi);
622 ts_sat.emplace_back(rhit.T()-
tns_min);
623 terrs_sat.emplace_back(st);
627 pos_unsat.emplace_back(
pos-Pi);
628 ts_unsat.emplace_back(rhit.T()-
tns_min);
629 terrs_unsat.emplace_back(st);
632 if (rhit.T() <
t_min ) {
639 if (rhit.T() >
t_max ) {
648 double ADC_sq_sum = 0;
649 for (
unsigned i = 0;
i != track.
NCell(); ++
i) {
651 double tns = chit->
TNS();
653 ADC_sq_sum+=chit->
ADC()*chit->
ADC();
657 adc_Center = adc_Center*(1/double(
ADC_tot));
658 GeV_Center = GeV_Center*(1/double(
GeV_tot));
659 geo_Center = geo_Center*(1/double(ts_sat.size()+ts_unsat.size()));
661 adc_W = (adc_Center-geo_Center).
Mag();
662 GeV_W = (GeV_Center - geo_Center).
Mag();
679 if (pos_unsat.size()<2) {
716 double max_time_gap = 0;
717 double time_cross =
hitxt.back()-
hitxt.front()+0.00000001;
720 if (dt > max_time_gap) max_time_gap =
dt;
725 time_cross =
hityt.back()-
hityt.front()+0.00000001;
728 if (dt > max_time_gap) max_time_gap =
dt;
746 const auto n = x_val.size();
747 const auto x = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/
n;
748 const auto y = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/
n;
749 const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/
n;
750 const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/
n;
751 const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/
n;
753 const auto b = (xy -
x*
y)/(
xx -
x*
x);
754 const auto a =
y -
b*
x;
755 const auto r = (xy - x*
y)/
sqrt((
xx - x*x)*(yy -
y*
y));
763 if (
fabs(x-x0) < 0.1 &&
fabs(y-y0)<0.1)
return 0;
764 double cos2 = (x-x0 + dydx*(y-
y0))*(x-x0 + dydx*(y-
y0))/
dist2(x-x0,y-y0)/
dist2(1,dydx);
765 double docasqr =
dist2(x-x0,y-y0)*(1-cos2);
771 std::vector<double> & Ts,
772 std::vector<double> & Terr,
784 std::vector<double> Vs;
786 for (
int i = 0;
i!=
m; ++
i ) {
787 sumL += Positions[
i].Mag();
789 sumLT += Positions[
i].Mag() * Ts[
i];
790 sumLL += Positions[
i].Mag() * Positions[
i].Mag();
793 Tmean = sumT/double(m);
794 double v = (sumL*sumL-m*sumLL)/(sumT*sumL-m*sumLT);
805 for (
int i = 0;
i!=
m; ++
i ) {
809 double ti =
genR->Gaus(Ts[
i], Terr[i]);
810 double li =
sqrt(xi*xi+yi*yi+zi*zi);
817 v = (sumL*sumL-m*sumLL)/(sumT*sumL-m*sumLT);
826 for (
unsigned j = 0;
j!=Vs.size(); ++
j)
828 v_mean /= double(Vs.size());
830 for (
unsigned j = 0;
j!=Vs.size(); ++
j)
831 v_err += (Vs[
j]-v_mean)*(Vs[
j]-v_mean);
833 v_err =
sqrt(v_err)/double(Vs.size());
835 if (v_mean < 0.003 && v_mean > 0)
838 else if (v_mean > -0.003 && v_mean < 0)
841 double t0 = Tmean-sumL/(m+0.)/v_mean;
845 for (
int i = 0;
i!=
m; ++
i ) {
846 Trms += (Ts[
i]-Tmean)*(Ts[
i]-Tmean);
847 dTsqr+= (Ts[
i]-t0-Positions[
i].Mag()/v_mean)* (Ts[
i]-t0-Positions[
i].
Mag()/v_mean);
849 Trms =
sqrt(Trms/(m+0.));
850 dTsqr =
sqrt(dTsqr/(m+0.));
T max(const caf::Proxy< T > &a, T b)
back track the reconstruction to the simulation
double GetTimeRes(rb::CellHit const &cellhit)
void analyze(const art::Event &evt)
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
double Py(const int i=0) const
bool isfinite(const stan::math::var &v)
fvar< T > fabs(const fvar< T > &x)
std::vector< int > hityPlane
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
std::vector< double > hitxx
void LocalToWorld(const double *local, double *world) const
unsigned short Plane() const
FmmTrackerValidation(fhicl::ParameterSet const &pset)
const CellGeo * Cell(int icell) const
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
list_type::const_iterator const_iterator
std::vector< int > hitxPlane
Vertical planes which measure X.
double Px(const int i=0) const
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
void LinFit(const std::vector< double > &x, const std::vector< double > &y, double *fitpar)
std::string fTrackerInput
A rb::Prong with full reconstructed trajectory.
std::vector< double > hityt
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
std::vector< int > hitxCell
Horizontal planes which measure Y.
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
unsigned short Cell() const
std::vector< double > hitxe
double P(const int i=0) const
Collect Geo headers and supply basic geometry functions.
double vfit(std::vector< TVector3 > &Positions, std::vector< double > &Ts, std::vector< double > &Terr, double &v_err, double &Tmean, double &Trms, double &dTsqr)
std::vector< double > hitxz
EDAnalyzer(Table< Config > const &config)
std::vector< double > hitye
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
double Vx(const int i=0) const
T * make(ARGS...args) const
double dist2(double dx, double dy)
int16_t ADC(uint32_t i) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
std::vector< double > hityz
double Pz(const int i=0) const
double doca2(double x, double y, double x0, double y0, double dydx)
virtual ~FmmTrackerValidation()
double Vz(const int i=0) const
std::vector< double > hitxt
std::vector< int > hityCell
T min(const caf::Proxy< T > &a, T b)
bool is_trigger_by_epoch1_fmmtrigger(art::Handle< std::vector< rb::Cluster > > slices, int &_penetrated, double &_nSatHits, double &_meanADC)
Encapsulate the cell geometry.
std::vector< double > hityy
double Vy(const int i=0) const
Encapsulate the geometry of one entire detector (near, far, ndos)