21 #include "Utilities/func/MathUtil.h" 30 #include "NovaDAQConventions/DAQConventions.h" 100 std::map<std::string, double>
t_;
101 std::map<std::string, double>
t2_;
107 void NumberOfCellsPerLength(
rb::Cluster const &
hits,
double tracklength_of_xview,
double tracklength_of_yview,
double &number_of_cells_per_length_xview,
double &number_of_cells_per_length_yview)
const;
108 float Distance(
unsigned const &xcell1,
unsigned const &ycell1,
unsigned const &plane1,
unsigned const &xcell2,
unsigned const &ycell2 ,
unsigned const &plane2)
const;
111 double WeightedCenterCut(
rb::Cluster const &
hits,
float PX_min,
float PX_max,
float PY_min,
float PY_max,
float CX_min,
float CX_max,
float CY_min,
float CY_max,
double &weighted_off_center_xx,
double &weighted_off_center_xz,
double &weighted_off_center_yy,
double &weighted_off_center_yz)
const;
114 void LinFit(
const std::vector<double>&
x,
const std::vector<double>&
y,
double *fitpar);
155 if (p->
PdgCode()!=42 )
continue;
158 unsigned ntruehits_Tot = flshits.size();
159 if ( ntruehits_Tot==0 )
continue;
161 double momentum2 = (p->
P())*(p->
P());
162 double mass2 = (p->
Mass())*(p->
Mass());
164 ti_MC = flshits[0].GetEntryT();
165 tf_MC = flshits[0].GetExitT();
166 size_t planeID_i = flshits[0].GetPlaneID();
169 for (
unsigned h =1;
h!= flshits.size(); ++
h) {
170 double ti = flshits[
h].GetEntryT();
171 double tf = flshits[
h].GetExitT();
195 for (
unsigned i = 0;
i!= slices->size(); ++
i) {
215 bool PassThrough =
false;
216 bool Intrusion =
false;
224 std::vector<double> x_hit;
225 std::vector<double> y_hit;
226 std::vector<double> zy_hit;
227 std::vector<double> zx_hit;
229 for (
auto const hit_ptr: hits.
AllCells()) {
232 sq_sum_ADC+=hit.
ADC()*hit.
ADC();
239 x_hit.push_back(hit.
Cell());
240 zx_hit.push_back(hit.
Plane());
258 else if (!PassThrough) {
260 else if (hit.
Cell()>
_xMax-
_deltx && entry_id != 2) PassThrough =
true;
267 y_hit.push_back(hit.
Cell());
268 zy_hit.push_back(hit.
Plane());
286 else if (!PassThrough) {
288 else if (hit.
Cell()>
_yMax-
_delty && entry_id != 4) PassThrough =
true;
302 else if (!PassThrough) {
310 double preselection=0;
330 double number_of_cells_per_length_xview=0;
331 double number_of_cells_per_length_yview=0;
332 NumberOfCellsPerLength(hits,tracklength_of_xview,tracklength_of_yview,number_of_cells_per_length_xview,number_of_cells_per_length_yview);
334 LinFit(x_hit, zx_hit, fitpar);
335 double r2x = fitpar[2];
336 LinFit(y_hit, zy_hit, fitpar);
337 double r2y = fitpar[2];
345 stdevADC =
std::sqrt( sq_sum_ADC*1.0/(nx+ny) - sumADC*1.0/(nx+ny)*sumADC*1.0/(nx+ny));
351 double weighted_off_center_xx=0, weighted_off_center_xz=0, weighted_off_center_yy=0, weighted_off_center_yz=0;
352 double stdev_cells_per_plane_xview=0, stdev_cells_per_plane_yview=0;
356 tset(
"slice_number",
i);
359 tset(
"number_of_hits_per_length",hits.
NCell()/tracklength);
360 tset(
"number_of_cells_per_length_xview",number_of_cells_per_length_xview);
361 tset(
"number_of_cells_per_length_yview",number_of_cells_per_length_yview);
362 tset(
"track_length",tracklength);
364 tset(
"stdev_cells_per_plane_xview",stdev_cells_per_plane_xview);
365 tset(
"stdev_cells_per_plane_yview",stdev_cells_per_plane_yview);
366 tset(
"weighted_off_center",
WeightedCenterCut(hits,
PX_min,
PX_max,
PY_min,
PY_max,
CX_min,
CX_max,
CY_min,
CY_max, weighted_off_center_xx, weighted_off_center_xz, weighted_off_center_yy, weighted_off_center_yz ));
367 tset(
"weighted_off_center_xx",weighted_off_center_xx);
368 tset(
"weighted_off_center_xz",weighted_off_center_xz);
369 tset(
"weighted_off_center_yy",weighted_off_center_yy);
370 tset(
"weighted_off_center_yz",weighted_off_center_yz);
371 tset(
"meanADC",sumADC*1.0/(nx+ny));
372 tset(
"stdevADC",stdevADC);
382 tset(
"preselection",preselection);
402 std::vector<std::string> branch_names =
403 {
"run_number",
"subrun_number",
"event_number",
"slice_number",
"number_of_hits_per_length",
"number_of_cells_per_length_xview",
"number_of_cells_per_length_yview",
"track_length",
"stdev_cells_per_plane_xview",
"stdev_cells_per_plane_yview",
"weighted_off_center",
"weighted_off_center_xx",
"weighted_off_center_xz",
"weighted_off_center_yy",
"weighted_off_center_yz",
"percent_of_hits_in_overlap_planes",
"number_of_surface_hits",
"betaMC",
"xhitsMC",
"yhitsMC",
"trueHits",
"preselection",
"meanADC",
"stdevADC",
"nxhits",
"nyhits",
"deltaTDC",
"deltaP",
"nSatHits",
"r2x",
"r2y",
"dPX",
"dPY",
"dCX",
"dCY"};
404 for (
auto const&
name : branch_names)
407 std::vector<std::string> branch_names2 =
408 {
"run_number",
"subrun_number",
"event_number",
"betaMC"};
409 for (
auto const&
name : branch_names2)
415 auto it =
t_.find(branch_name);
418 std::cerr <<
"\n\t Branch " << branch_name <<
" does not exist! \n" 428 auto it =
t2_.find(branch_name);
431 std::cerr <<
"\n\t Branch " << branch_name <<
" does not exist! \n" 440 std::vector<std::vector<int>> x_hits_plane(
PX_max-
PX_min+1,std::vector<int>(0));
441 std::vector<std::vector<int>> y_hits_plane(
PY_max-
PY_min+1,std::vector<int>(0));
442 for (
auto const hit_ptr: hits.
AllCells()) {
456 for (
unsigned i=0;
i!=x_hits_plane.size();
i++){
457 std::sort(x_hits_plane[
i].begin(),x_hits_plane[
i].
end());
458 std::vector<int>::iterator
it;
459 it = std::unique (x_hits_plane[
i].begin(),x_hits_plane[
i].
end());
461 x_cells_plane[
i]=x_hits_plane[
i].size();
463 for (
unsigned i=0;
i!=y_hits_plane.size();
i++){
464 std::sort(y_hits_plane[
i].begin(),y_hits_plane[
i].
end());
465 std::vector<int>::iterator
it;
466 it = std::unique (y_hits_plane[
i].begin(),y_hits_plane[
i].
end());
468 y_cells_plane[
i]=y_hits_plane[
i].size();
472 for (
unsigned i=0;
i!=x_cells_plane.size();
i++){
473 sum_of_cells+=x_cells_plane[
i];
475 number_of_cells_per_length_xview=sum_of_cells/tracklength_of_xview;
477 for (
unsigned i=0;
i!=y_cells_plane.size();
i++){
478 sum_of_cells+=y_cells_plane[
i];
480 number_of_cells_per_length_yview=sum_of_cells/tracklength_of_yview;
487 std::vector<std::vector<int>> x_hits_plane(
PX_max-
PX_min+1,std::vector<int>(0));
488 std::vector<std::vector<int>> y_hits_plane(
PY_max-
PY_min+1,std::vector<int>(0));
489 for (
auto const hit_ptr: hits.
AllCells()) {
507 for (
unsigned i=0;
i!=x_hits_plane.size();
i++){
508 std::sort(x_hits_plane[
i].begin(),x_hits_plane[
i].
end());
509 std::vector<int>::iterator
it;
510 it = std::unique (x_hits_plane[
i].begin(),x_hits_plane[
i].
end());
512 x_cells_plane[
i]=x_hits_plane[
i].size();
514 for (
unsigned i=0;
i!=y_hits_plane.size();
i++){
515 std::sort(y_hits_plane[
i].begin(),y_hits_plane[
i].
end());
516 std::vector<int>::iterator
it;
517 it = std::unique (y_hits_plane[
i].begin(),y_hits_plane[
i].
end());
519 y_cells_plane[
i]=y_hits_plane[
i].size();
521 sum = std::accumulate(x_cells_plane.begin(),x_cells_plane.end(),0.0);
522 mean = sum / x_cells_plane.size();
523 sq_sum = std::inner_product(x_cells_plane.begin(),x_cells_plane.end(),x_cells_plane.begin(),0.0);
524 stdev =
std::sqrt(sq_sum / x_cells_plane.size()-mean *
mean);
527 stdev_cells_per_plane_xview=stdev/
mean;
528 sum = std::accumulate(y_cells_plane.begin(),y_cells_plane.end(),0.0);
529 mean = sum / y_cells_plane.size();
530 sq_sum = std::inner_product(y_cells_plane.begin(),y_cells_plane.end(),y_cells_plane.begin(),0.0);
531 stdev =
std::sqrt(sq_sum / y_cells_plane.size()-mean *
mean);
532 stdev_cells_per_plane_yview=stdev/
mean;
537 float zcl::FastMMStudy::Distance(
unsigned const &xcell1,
unsigned const &ycell1,
unsigned const &plane1,
unsigned const &xcell2,
unsigned const &ycell2 ,
unsigned const &plane2)
const {
538 return std::sqrt((plane1-plane2)*(plane1-plane2)*2.82072+(xcell1-xcell2)*(xcell1-xcell2)+(ycell1-ycell2)*(ycell1-ycell2));
584 int number_of_surface_hits=0;
585 for(
unsigned i=0;
i!=hits.
NCell();++
i){
588 return number_of_surface_hits;
590 double zcl::FastMMStudy::WeightedCenterCut(
rb::Cluster const &
hits,
float PX_min,
float PX_max,
float PY_min,
float PY_max,
float CX_min,
float CX_max,
float CY_min,
float CY_max,
double &weighted_off_center_xx,
double &weighted_off_center_xz,
double &weighted_off_center_yy,
double &weighted_off_center_yz )
const{
591 float Weighted_PX_Center=0.5,Weighted_CX_Center=0.5,Weighted_PY_Center=0.5,Weighted_CY_Center=0.5,SumOfXHitsADC=0,SumOfYHitsADC=0;
592 for (
auto const hit_ptr: hits.
AllCells()) {
595 SumOfXHitsADC+=hit.
ADC();
600 SumOfYHitsADC+=hit.
ADC();
605 if(PX_max!=PX_min) Weighted_PX_Center=Weighted_PX_Center/SumOfXHitsADC/(PX_max-
PX_min);
606 if(PY_max!=PY_min) Weighted_PY_Center=Weighted_PY_Center/SumOfYHitsADC/(PY_max-
PY_min);
607 if(CX_max!=CX_min) Weighted_CX_Center=Weighted_CX_Center/SumOfXHitsADC/(CX_max-
CX_min);
608 if(CY_max!=CY_min) Weighted_CY_Center=Weighted_CY_Center/SumOfYHitsADC/(CY_max-
CY_min);
612 weighted_off_center_xx=
std::abs(Weighted_CX_Center-0.5);
613 weighted_off_center_xz=
std::abs(Weighted_PX_Center-0.5);
614 weighted_off_center_yy=
std::abs(Weighted_CY_Center-0.5);
615 weighted_off_center_yz=
std::abs(Weighted_PY_Center-0.5);
622 double number_hits_in_overlap_planes=0;
623 for (
auto const hit_ptr: hits.
AllCells()) {
625 if(hit.
Plane()<=max && hit.
Plane()>=
min) number_hits_in_overlap_planes++;
628 return number_hits_in_overlap_planes/hits.
NCell();
634 const auto n = x_val.size();
635 const auto x = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/
n;
636 const auto y = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/
n;
637 const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/
n;
638 const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/
n;
639 const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/
n;
641 const auto b = (xy -
x*
y)/(
xx -
x*
x);
642 const auto a =
y -
b*
x;
643 const auto r = (xy - x*
y)/
sqrt((
xx - x*x)*(yy -
y*
y));
std::map< std::string, double > t_
void StdevCellsPerPlane(rb::Cluster const &hits, double &stdev_cells_per_plane_xview, double &stdev_cells_per_plane_yview) const
T max(const caf::Proxy< T > &a, T b)
SubRunNumber_t subRun() const
back track the reconstruction to the simulation
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
double branch_default_value_
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
int32_t TDC() const
The time of the last baseline sample.
void analyze(art::Event const &e) override
FastMMStudy & operator=(FastMMStudy const &)=delete
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
std::string _sliceinstance
void LinFit(const std::vector< double > &x, const std::vector< double > &y, double *fitpar)
double _stdevcellsperplane
unsigned short Plane() 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
Vertical planes which measure X.
void tset(std::string const &branch_name, double const &value)
A collection of associated CellHits.
double WeightedCenterCut(rb::Cluster const &hits, float PX_min, float PX_max, float PY_min, float PY_max, float CX_min, float CX_max, float CY_min, float CY_max, double &weighted_off_center_xx, double &weighted_off_center_xz, double &weighted_off_center_yy, double &weighted_off_center_yz) const
double NumberOfHitsInOverlapPlanesCut(rb::Cluster const &hits, float PX_min, float PX_max, float PY_min, float PY_max) const
DEFINE_ART_MODULE(TestTMapFile)
unsigned distance(const T &t1, const T &t2)
Horizontal planes which measure Y.
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
bool SurfAssign(rb::CellHit const &hit) const
unsigned short Cell() const
art::ServiceHandle< art::TFileService > tfs_
const XML_Char int const XML_Char * value
int _MinNumberSurfaceHits
double P(const int i=0) const
Collect Geo headers and supply basic geometry functions.
EventNumber_t event() const
EDAnalyzer(Table< Config > const &config)
static float min(const float a, const float b, const float c)
void NumberOfCellsPerLength(rb::Cluster const &hits, double tracklength_of_xview, double tracklength_of_yview, double &number_of_cells_per_length_xview, double &number_of_cells_per_length_yview) const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
int NumberOfSurfaceHits(rb::Cluster const &hits) const
A rawdata::RawDigit with channel information decoded.
T * make(ARGS...args) const
int16_t ADC(uint32_t i) const
void tset2(std::string const &branch_name, double const &value)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
int GetSurfaceId(rb::CellHit const &hit) const
assert(nhit_max >=nhit_nbins)
T min(const caf::Proxy< T > &a, T b)
float Distance(unsigned const &xcell1, unsigned const &ycell1, unsigned const &plane1, unsigned const &xcell2, unsigned const &ycell2, unsigned const &plane2) const
std::map< std::string, double > t2_
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
FastMMStudy(fhicl::ParameterSet const &p)
Encapsulate the geometry of one entire detector (near, far, ndos)