Public Member Functions | Public Attributes | Protected Member Functions | Static Protected Member Functions | Protected Attributes | List of all members
cluster::ClusterParamsAlg Class Reference

#include <ClusterParamsAlg.h>

Public Member Functions

 ClusterParamsAlg ()
 
 ClusterParamsAlg (const std::vector< util::PxHit > &)
 
void Initialize ()
 
void SetMinNHits (size_t nhit)
 
size_t MinNHits () const
 
int SetHits (const std::vector< util::PxHit > &)
 
void SetRefineDirectionQMin (double qmin)
 
void SetVerbose (bool yes=true)
 
template<typename Stream >
void TimeReport (Stream &stream) const
 
void GetFANNVector (std::vector< float > &data)
 
void PrintFANNVector ()
 
void FillParams (util::GeometryUtilities const &gser, bool override_DoGetAverages=false, bool override_DoGetRoughAxis=false, bool override_DoGetProfileInfo=false, bool override_DoRefineStartPointsAndDirection=false, bool override_DoGetFinalSlope=false, bool override_DoTrackShowerSep=false, bool override_DoEndCharge=false)
 
const cluster_paramsGetParams () const
 
void GetAverages (bool override=false)
 
void GetRoughAxis (bool override=false)
 
void GetProfileInfo (util::GeometryUtilities const &gser, bool override=false)
 
void RefineStartPoints (util::GeometryUtilities const &gser)
 
void GetFinalSlope (util::GeometryUtilities const &gser, bool override=false)
 
void GetEndCharges (util::GeometryUtilities const &gser, bool override_=false)
 
void RefineDirection (bool override=false)
 
void RefineStartPointAndDirection (util::GeometryUtilities const &gser, bool override=false)
 
void TrackShowerSeparation (bool override=false)
 
void setNeuralNetPath (std::string s)
 
void FillPolygon (util::GeometryUtilities const &gser)
 
void GetOpeningAngle ()
 
const util::PxPointRoughStartPoint ()
 
const util::PxPointRoughEndPoint ()
 
double RoughSlope ()
 
double RoughIntercept ()
 
double StartCharge (util::GeometryUtilities const &gser, float length=1., unsigned int nbins=10)
 Returns the expected charge at the beginning of the cluster. More...
 
double EndCharge (util::GeometryUtilities const &gser, float length=1., unsigned int nbins=10)
 Returns the expected charge at the end of the cluster. More...
 
float MultipleHitWires ()
 Returns the number of multiple hits per wire. More...
 
float MultipleHitDensity (util::GeometryUtilities const &gser)
 Returns the number of multiple hits per wire. More...
 
void EnableFANN ()
 
void DisableFANN ()
 
size_t GetNHits () const
 
const std::vector< util::PxHit > & GetHitVector () const
 
int Plane () const
 
void SetPlane (int p)
 

Public Attributes

cluster::cluster_params fParams
 
std::string fNeuralNetPath
 
std::vector< std::stringfTimeRecord_ProcName
 
std::vector< double > fTimeRecord_ProcTime
 

Protected Member Functions

double IntegrateFitCharge (util::GeometryUtilities const &gser, double from_length, double to_length, unsigned int fit_first_bin, unsigned int fit_end_bin)
 Integrates the charge between two positions in the cluster axis. More...
 

Static Protected Member Functions

static double LinearIntegral (double m, double q, double x1, double x2)
 Returns the integral of f(x) = mx + q defined in [x1, x2]. More...
 

Protected Attributes

size_t fMinNHits
 Cut value for # hits: below this value clusters are not evaluated. More...
 
std::vector< util::PxHitfHitVector
 
bool verbose
 
std::vector< double > fChargeCutoffThreshold
 
int fPlane
 
double fQMinRefDir
 
std::vector< double > fChargeProfile
 
std::vector< double > fCoarseChargeProfile
 
std::vector< double > fChargeProfileNew
 
int fCoarseNbins
 
int fProfileNbins
 
int fProfileMaximumBin
 
double fProfileIntegralForward
 
double fProfileIntegralBackward
 
double fProjectedLength
 
double fBeginIntercept
 
double fEndIntercept
 
double fInterHigh_side
 
double fInterLow_side
 
bool fFinishedGetAverages
 
bool fFinishedGetRoughAxis
 
bool fFinishedGetProfileInfo
 
bool fFinishedRefineStartPoints
 
bool fFinishedRefineDirection
 
bool fFinishedGetFinalSlope
 
bool fFinishedRefineStartPointAndDirection
 
bool fFinishedTrackShowerSep
 
bool fFinishedGetEndCharges
 
double fRough2DSlope
 
double fRough2DIntercept
 
util::PxPoint fRoughBeginPoint
 
util::PxPoint fRoughEndPoint
 
bool enableFANN
 

Detailed Description

Definition at line 24 of file ClusterParamsAlg.h.

Constructor & Destructor Documentation

cluster::ClusterParamsAlg::ClusterParamsAlg ( )

Definition at line 34 of file ClusterParamsAlg.cxx.

35  {
36  fMinNHits = 10;
37  enableFANN = false;
38  verbose = true;
39  Initialize();
40  }
size_t fMinNHits
Cut value for # hits: below this value clusters are not evaluated.
cluster::ClusterParamsAlg::ClusterParamsAlg ( const std::vector< util::PxHit > &  inhitlist)

Definition at line 42 of file ClusterParamsAlg.cxx.

43  {
44  fMinNHits = 10;
45  enableFANN = false;
46  verbose = true;
47  SetHits(inhitlist);
48  }
size_t fMinNHits
Cut value for # hits: below this value clusters are not evaluated.
int SetHits(const std::vector< util::PxHit > &)

Member Function Documentation

void cluster::ClusterParamsAlg::DisableFANN ( )
inline

Definition at line 268 of file ClusterParamsAlg.h.

269  {
270  enableFANN = false;
271  }
void cluster::ClusterParamsAlg::EnableFANN ( )

Definition at line 180 of file ClusterParamsAlg.cxx.

181  {
182  enableFANN = true;
183  }
double cluster::ClusterParamsAlg::EndCharge ( util::GeometryUtilities const &  gser,
float  length = 1.,
unsigned int  nbins = 10 
)

Returns the expected charge at the end of the cluster.

Parameters
nbinsuse at least this number of charge bins from charge profile
lengthspace before the end of cluster where to collect charge, in cm the expected charge at the end of the cluster
See also
StartCharge(), IntegrateFitCharge()

This method returns the charge under the last length cm of the cluster. See StartCharge() for a detailed explanation. For even more details, see IntegrateFitCharge().

Definition at line 1487 of file ClusterParamsAlg.cxx.

1490  {
1491  switch (fHitVector.size()) {
1492  case 0: return 0.;
1493  case 1:
1494  return fHitVector.back().charge;
1495  // the "default" is the rest of the function
1496  } // switch
1497 
1498  // need the available number of bins and the axis length
1499  GetProfileInfo(gser);
1500  const unsigned int MaxBins = fChargeProfile.size();
1501 
1502  // this is the range of the fit:
1503  const unsigned int fit_first_bin = MaxBins > nbins ? MaxBins - nbins : 0,
1504  fit_last_bin = MaxBins;
1505 
1506  // now determine the integration range, in bin units;
1507  // get to the end, and go length backward;
1508  // note that length can be pathologic (0, negative...); not our problem!
1509  const double from = fProjectedLength - length, to = fProjectedLength;
1510 
1511  return IntegrateFitCharge(gser, from, to, fit_first_bin, fit_last_bin);
1512  } // ClusterParamsAlg::EndCharge()
double IntegrateFitCharge(util::GeometryUtilities const &gser, double from_length, double to_length, unsigned int fit_first_bin, unsigned int fit_end_bin)
Integrates the charge between two positions in the cluster axis.
void GetProfileInfo(util::GeometryUtilities const &gser, bool override=false)
std::vector< util::PxHit > fHitVector
std::vector< double > fChargeProfile
void cluster::ClusterParamsAlg::FillParams ( util::GeometryUtilities const &  gser,
bool  override_DoGetAverages = false,
bool  override_DoGetRoughAxis = false,
bool  override_DoGetProfileInfo = false,
bool  override_DoRefineStartPointsAndDirection = false,
bool  override_DoGetFinalSlope = false,
bool  override_DoTrackShowerSep = false,
bool  override_DoEndCharge = false 
)

Runs all the functions which calculate cluster params and stashes the results in the private ClusterParams struct.

Parameters
override_DoGetAveragesforce re-execution of GetAverages()
override_DoGetRoughAxisforce re-execution of GetRoughAxis()
override_DoGetProfileInfoforce re-execution of GetProfileInfo()
override_DoRefineStartPointsforce re-execution of RefineStartPoints()
override_DoGetFinalSlopeforce re-execution of GetFinalSlope()
override_DoEndChargeforce re-execution of GetEndCharges()

Definition at line 186 of file ClusterParamsAlg.cxx.

194  {
195  GetAverages(override_DoGetAverages);
196  GetRoughAxis(override_DoGetRoughAxis);
197  GetProfileInfo(gser, override_DoGetProfileInfo);
198  RefineStartPointAndDirection(gser, override_DoStartPointsAndDirection);
199  GetFinalSlope(gser, override_DoGetFinalSlope);
200  GetEndCharges(gser, override_DoEndCharge);
201  TrackShowerSeparation(override_DoTrackShowerSep);
202  }
void GetRoughAxis(bool override=false)
void GetFinalSlope(util::GeometryUtilities const &gser, bool override=false)
void GetProfileInfo(util::GeometryUtilities const &gser, bool override=false)
void GetEndCharges(util::GeometryUtilities const &gser, bool override_=false)
void TrackShowerSeparation(bool override=false)
void GetAverages(bool override=false)
void RefineStartPointAndDirection(util::GeometryUtilities const &gser, bool override=false)
void cluster::ClusterParamsAlg::FillPolygon ( util::GeometryUtilities const &  gser)

Definition at line 1270 of file ClusterParamsAlg.cxx.

1271  {
1272  TStopwatch localWatch;
1273  localWatch.Start();
1274 
1275  if (fHitVector.size()) {
1276  std::vector<const util::PxHit*> container_polygon;
1277  gser.SelectPolygonHitList(fHitVector, container_polygon);
1278  //now making Polygon Object
1279  std::pair<float, float> tmpvertex;
1280  //make Polygon Object as in mac/PolyOverlap.cc
1281  std::vector<std::pair<float, float>> vertices;
1282  for (unsigned int i = 0; i < container_polygon.size(); i++) {
1283  tmpvertex = std::make_pair(container_polygon.at(i)->w, container_polygon.at(i)->t);
1284  vertices.push_back(tmpvertex);
1285  }
1286  fParams.PolyObject = Polygon2D(vertices);
1287  }
1288 
1289  fTimeRecord_ProcName.push_back("FillPolygon");
1290  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1291  }
std::vector< std::string > fTimeRecord_ProcName
Polygon2D PolyObject
Polygon Object...see Polygon2D.hh.
Definition: ClusterParams.h:22
std::vector< util::PxHit > fHitVector
cluster::cluster_params fParams
std::vector< double > fTimeRecord_ProcTime
void cluster::ClusterParamsAlg::GetAverages ( bool  override = false)

Calculates the following variables: mean_charge mean_x mean_y charge_wgt_x charge_wgt_y eigenvalue_principal eigenvalue_secondary multi_hit_wires N_Wires

Parameters
overrideforce recalculation of variables

Definition at line 205 of file ClusterParamsAlg.cxx.

206  {
207  if (!override) { //Override being set, we skip all this logic.
208  //OK, no override. Stop if we're already finshed.
209  if (fFinishedGetAverages) return;
210  }
211 
212  TStopwatch localWatch;
213  localWatch.Start();
214 
215  TPrincipal fPrincipal(2, "D");
216 
217  fParams.N_Hits = fHitVector.size();
218 
219  std::map<double, int> wireMap;
220 
221  lar::util::StatCollector<double> charge, sumADC;
222 
223  int uniquewires = 0;
224  int multi_hit_wires = 0;
225  for (auto& hit : fHitVector) {
226  double data[2];
227  data[0] = hit.w;
228  data[1] = hit.t;
229  fPrincipal.AddRow(data);
230  fParams.charge_wgt_x += hit.w * hit.charge;
231  fParams.charge_wgt_y += hit.t * hit.charge;
232  charge.add(hit.charge);
233  sumADC.add(hit.sumADC);
234 
235  if (wireMap[hit.w] == 0) { uniquewires++; }
236  if (wireMap[hit.w] == 1) { multi_hit_wires++; }
237  wireMap[hit.w]++;
238  }
239 
240  fParams.sum_charge = charge.Sum();
241  fParams.mean_charge = charge.Average();
242  fParams.rms_charge = charge.RMS();
243 
244  fParams.sum_ADC = sumADC.Sum();
245  fParams.mean_ADC = sumADC.Average();
246  fParams.rms_ADC = sumADC.RMS();
247 
248  fParams.N_Wires = uniquewires;
249  fParams.multi_hit_wires = multi_hit_wires;
250 
251  if (fPrincipal.GetMeanValues()->GetNrows() < 2) { throw cluster::CRUException(); }
252 
253  fParams.mean_x = (*fPrincipal.GetMeanValues())[0];
254  fParams.mean_y = (*fPrincipal.GetMeanValues())[1];
256 
257  if (fParams.sum_charge != 0.) {
260  }
261  else { // "SNAFU"; use the mean
264  }
265 
266  fPrincipal.MakePrincipals();
267 
268  fParams.eigenvalue_principal = (*fPrincipal.GetEigenValues())[0];
269  fParams.eigenvalue_secondary = (*fPrincipal.GetEigenValues())[1];
270 
271  fFinishedGetAverages = true;
272 
273  fTimeRecord_ProcName.push_back("GetAverages");
274  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
275  }
double rms_ADC
RMS (standard deviation of sample) of ADC counts of hits in ADC.
Definition: ClusterParams.h:32
std::vector< std::string > fTimeRecord_ProcName
double charge_wgt_y
Mean of hits along y, charge weighted.
Definition: ClusterParams.h:38
double mean_x
Mean of hits along x, peaks only.
Definition: ClusterParams.h:33
double rms_charge
RMS (standard deviation of sample) of charge of hits in ADC.
Definition: ClusterParams.h:29
std::vector< util::PxHit > fHitVector
double eigenvalue_principal
the principal eigenvalue from PCA
Definition: ClusterParams.h:47
double mean_ADC
Mean (average) of ADC counts of hits, in ADC.
Definition: ClusterParams.h:31
double mean_y
Mean of hits along y, peaks only.
Definition: ClusterParams.h:34
Weight_t RMS() const
Returns the root mean square.
Weight_t Average() const
Returns the value average.
double eigenvalue_secondary
the secondary eigenvalue from PCA
Definition: ClusterParams.h:48
double sum_charge
Sum charge of hits in ADC.
Definition: ClusterParams.h:27
cluster::cluster_params fParams
Weight_t Sum() const
Returns the weighted sum of the values.
std::vector< double > fTimeRecord_ProcTime
Detector simulation of raw signals on wires.
double sum_ADC
Sum charge of ADC counts of hits, in ADC.
Definition: ClusterParams.h:30
double mean_charge
Mean (average) charge of hits in ADC.
Definition: ClusterParams.h:28
double charge_wgt_x
Mean of hits along x, charge weighted.
Definition: ClusterParams.h:37
Collects statistics on a single quantity (weighted)
void add(Data_t value, Weight_t weight=Weight_t(1.0))
Adds one entry with specified value and weight.
void cluster::ClusterParamsAlg::GetEndCharges ( util::GeometryUtilities const &  gser,
bool  override_ = false 
)

Calculates the following variables: start_charge end_charge

Parameters
override_force recompute the variables
See also
StartCharge(), EndCharge()

Definition at line 1380 of file ClusterParamsAlg.cxx.

1381  {
1382  if (!override_) { //Override being set, we skip all this logic.
1383  //OK, no override. Stop if we're already finshed.
1384  if (fFinishedGetEndCharges) return;
1385  }
1386 
1387  TStopwatch localWatch;
1388  localWatch.Start();
1389 
1391  fParams.end_charge = EndCharge(gser);
1392 
1393  fFinishedGetEndCharges = true;
1394 
1395  fTimeRecord_ProcName.push_back("GetEndCharges");
1396  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1397  } // ClusterParamsAlg::GetEndCharges()
std::vector< std::string > fTimeRecord_ProcName
double start_charge
Charge at the start of the cluster.
Definition: ClusterParams.h:45
cluster::cluster_params fParams
std::vector< double > fTimeRecord_ProcTime
double EndCharge(util::GeometryUtilities const &gser, float length=1., unsigned int nbins=10)
Returns the expected charge at the end of the cluster.
double StartCharge(util::GeometryUtilities const &gser, float length=1., unsigned int nbins=10)
Returns the expected charge at the beginning of the cluster.
double end_charge
Charge at the (other) end of the cluster.
Definition: ClusterParams.h:46
void cluster::ClusterParamsAlg::GetFANNVector ( std::vector< float > &  data)

This function returns a feature vector suitable for a neural net This function uses the data from cluster_params but packages it up in a different way, and so is inappropriate to include in clusterParams.hh. That's why it's here.

Parameters
datatakes a reference to a vector< float>

Definition at line 89 of file ClusterParamsAlg.cxx.

90  {
91  unsigned int length = 13;
92  if (data.size() != length) data.resize(length);
93  data[0] = fParams.opening_angle / PI;
95  data[2] = fParams.closing_angle / PI;
97  data[4] = fParams.eigenvalue_principal;
98  data[5] = fParams.eigenvalue_secondary;
99  data[6] = fParams.width / fParams.length;
102  data[9] = fParams.N_Hits_HC / fParams.N_Hits;
103  data[10] = fParams.modified_hit_density;
104  data[11] = fParams.RMS_charge / fParams.mean_charge;
105  data[12] = log(fParams.sum_charge / fParams.length);
106  return;
107  }
double closing_angle_charge_wgt
Same for charge_wgt.
Definition: ClusterParams.h:44
double eigenvalue_principal
the principal eigenvalue from PCA
Definition: ClusterParams.h:47
constexpr double PI
double eigenvalue_secondary
the secondary eigenvalue from PCA
Definition: ClusterParams.h:48
double sum_charge
Sum charge of hits in ADC.
Definition: ClusterParams.h:27
cluster::cluster_params fParams
double mean_charge
Mean (average) charge of hits in ADC.
Definition: ClusterParams.h:28
double opening_angle_charge_wgt
Same for charge_wgt.
Definition: ClusterParams.h:42
double opening_angle
Width of angular distubtion wrt vertx.
Definition: ClusterParams.h:41
double closing_angle
Width of angular distubtion wrt endpoint.
Definition: ClusterParams.h:43
void cluster::ClusterParamsAlg::GetFinalSlope ( util::GeometryUtilities const &  gser,
bool  override = false 
)

Calculates the following variables: hit_density_1D hit_density_2D angle_2d direction

Parameters
override[description]

Calculates the following variables: angle_2d modified_hit_density

end testing

Definition at line 894 of file ClusterParamsAlg.cxx.

895  {
896  if (!override) { //Override being set, we skip all this logic.
897  //OK, no override. Stop if we're already finshed.
898  if (fFinishedGetFinalSlope) return;
899  //Try to run the previous function if not yet done.
901  }
902  else {
903  //Try to run the previous function if not yet done.
905  }
906 
907  TStopwatch localWatch;
908  localWatch.Start();
909  /**
910  * Calculates the following variables:
911  * angle_2d
912  * modified_hit_density
913  */
914 
915  const int NBINS = 720;
916  std::vector<int> fh_omega_single(NBINS, 0); //720,-180., 180.
917 
918  double current_maximum = 0;
919  double curr_max_bin = -1;
920 
921  for (auto hit : fHitVector) {
922 
923  double omx = gser.Get2Dangle((util::PxPoint*)&hit,
924  &fParams.start_point); // in rad and assuming cm/cm space.
925  int nbin = (omx + TMath::Pi()) * (NBINS - 1) / (2 * TMath::Pi());
926  if (nbin >= NBINS) nbin = NBINS - 1;
927  if (nbin < 0) nbin = 0;
928  fh_omega_single[nbin] += hit.charge;
929 
930  if (fh_omega_single[nbin] > current_maximum) {
931  current_maximum = fh_omega_single[nbin];
932  curr_max_bin = nbin;
933  }
934  }
935 
936  fParams.angle_2d = (curr_max_bin / 720 * (2 * TMath::Pi())) - TMath::Pi();
937  fParams.angle_2d *= 180 / PI;
938  if (verbose) std::cout << " Final 2D angle: " << fParams.angle_2d << " degrees " << std::endl;
939 
940  double mod_angle =
941  (fabs(fParams.angle_2d) <= 90) ?
942  fabs(fParams.angle_2d) :
943  180 - fabs(fParams.angle_2d); //want to transfer angle to radians and from 0 to 90.
944 
945  if (
946  cos(
947  mod_angle * PI /
948  180.)) { //values here based on fit of hit_density_1D vs. mod_angle defined as above (in ArgoNeuT).
949  if (mod_angle <= 75.)
951  fParams.hit_density_1D / (2.45 * cos(mod_angle * PI / 180.) + 0.2595);
952  else
954  fParams.hit_density_1D / (1.036 * mod_angle * PI / 180. - 0.2561);
955 
956  //calculate modified mean_charge:
957  fParams.modmeancharge = fParams.mean_charge / (1264 - 780 * cos(mod_angle * PI / 180.));
958  }
959  else
961 
962  /////////////////// testing
963  // do not use for now.
964  bool drawProfileHisto = false;
965  if (drawProfileHisto) {
966 
968  double corr_factor = 1;
969  if (
970  cos(
971  mod_angle * PI /
972  180.)) { //values here based on fit of hit_density_1D vs. mod_angle defined as above (in ArgoNeuT).
973 
974  if (mod_angle <= 75.)
975  corr_factor *= (2.45 * cos(mod_angle * PI / 180.) + 0.2595);
976  else
977  corr_factor *= (1.036 * mod_angle * PI / 180. - 0.2561);
978  }
979 
980  fProfileNbins =
981  (int)(fProfileNbins / 2 * corr_factor + 0.5); // +0.5 to round to nearest sensible value
982 
983  if (verbose) std::cout << " number of final profile bins " << fProfileNbins << std::endl;
984 
985  fChargeProfile.clear();
986  fCoarseChargeProfile.clear();
987  fChargeProfile.resize(fProfileNbins, 0);
989 
990  fChargeProfileNew.clear();
991  fChargeProfileNew.resize(fProfileNbins, 0);
992 
993  util::PxPoint BeginOnlinePoint;
994  gser.GetPointOnLine(fRough2DSlope, fRough2DIntercept, &fParams.start_point, BeginOnlinePoint);
995 
996  for (auto hit : fHitVector) {
997 
998  util::PxPoint OnlinePoint;
999  // get coordinates of point on axis.
1000  gser.GetPointOnLine(fRough2DSlope, &BeginOnlinePoint, &hit, OnlinePoint);
1001 
1002  double linedist = gser.Get2DDistance(&OnlinePoint, &BeginOnlinePoint);
1003  int fine_bin = (int)(linedist / fProjectedLength * (fProfileNbins - 1));
1004 
1005  if (fine_bin < fProfileNbins) //only fill if bin number is in range
1006  {
1007  fChargeProfileNew.at(fine_bin) += hit.charge;
1008  }
1009  }
1010 
1011  TH1F* charge_histo = new TH1F("charge dist", "charge dist", fProfileNbins, 0, fProfileNbins);
1012 
1013  TH1F* charge_diff = new TH1F("charge diff", "charge diff", fProfileNbins, 0, fProfileNbins);
1014 
1015  for (int ix = 0; ix < fProfileNbins - 1; ix++) {
1016  charge_histo->SetBinContent(ix, fChargeProfileNew[ix]);
1017 
1018  if (ix > 2 && ix < fProfileNbins - 3) {
1019  double diff =
1020  (1. / 12. * fChargeProfileNew[ix - 2] - 2. / 3. * fChargeProfileNew[ix - 1] +
1021  2. / 3. * fChargeProfileNew[ix + 1] - 1. / 12. * fChargeProfileNew[ix + 2]) /
1022  (double)fProfileNbins;
1023  charge_diff->SetBinContent(ix, diff);
1024  }
1025  }
1026 
1027  TCanvas* chCanv = new TCanvas("chCanv", "chCanv", 600, 600);
1028  chCanv->cd();
1029  charge_histo->SetLineColor(3);
1030  charge_histo->Draw("");
1031  charge_diff->SetLineColor(2);
1032  charge_diff->Draw("same");
1033  chCanv->Update();
1034  }
1035 
1036  /// end testing
1037 
1038  fTimeRecord_ProcName.push_back("GetFinalSlope");
1039  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1040 
1041  fFinishedGetFinalSlope = true;
1042  return;
1043  }
std::vector< double > fCoarseChargeProfile
std::vector< std::string > fTimeRecord_ProcName
std::vector< util::PxHit > fHitVector
constexpr double PI
void RefineStartPoints(util::GeometryUtilities const &gser)
cluster::cluster_params fParams
std::vector< double > fTimeRecord_ProcTime
util::PxPoint start_point
start point
Definition: ClusterParams.h:24
double angle_2d
Angle of axis in wire/hit view.
Definition: ClusterParams.h:40
Detector simulation of raw signals on wires.
std::vector< double > fChargeProfile
std::vector< double > fChargeProfileNew
double mean_charge
Mean (average) charge of hits in ADC.
Definition: ClusterParams.h:28
QTextStream & endl(QTextStream &s)
const std::vector<util::PxHit>& cluster::ClusterParamsAlg::GetHitVector ( ) const
inline

Definition at line 279 of file ClusterParamsAlg.h.

280  {
281  return fHitVector;
282  }
std::vector< util::PxHit > fHitVector
size_t cluster::ClusterParamsAlg::GetNHits ( ) const
inline

Definition at line 274 of file ClusterParamsAlg.h.

275  {
276  return fHitVector.size();
277  }
std::vector< util::PxHit > fHitVector
void cluster::ClusterParamsAlg::GetOpeningAngle ( )
const cluster_params& cluster::ClusterParamsAlg::GetParams ( ) const
inline

Definition at line 99 of file ClusterParamsAlg.h.

100  {
101  return fParams;
102  }
cluster::cluster_params fParams
void cluster::ClusterParamsAlg::GetProfileInfo ( util::GeometryUtilities const &  gser,
bool  override = false 
)

Calculates the following variables: opening_angle opening_angle_highcharge closing_angle closing_angle_highcharge offaxis_hits

Parameters
override[description]

Definition at line 349 of file ClusterParamsAlg.cxx.

350  {
351  if (!override) { //Override being set, we skip all this logic.
352  //OK, no override. Stop if we're already finshed.
353  if (fFinishedGetProfileInfo) return;
354  //Try to run the previous function if not yet done.
356  }
357  else {
359  }
360 
361  TStopwatch localWatch;
362  localWatch.Start();
363 
364  bool drawOrtHistos = false;
365 
366  //these variables need to be initialized to other values?
367  if (fRough2DSlope == -999.999 || fRough2DIntercept == -999.999) GetRoughAxis(true);
368 
369  //get slope of lines orthogonal to those found crossing the shower.
370  double inv_2d_slope = 0;
371  if (fabs(fRough2DSlope) > 0.001)
372  inv_2d_slope = -1. / fRough2DSlope;
373  else
374  inv_2d_slope = -999999.;
375 
376  double InterHigh = -999999;
377  double InterLow = 999999;
378  double WireHigh = -999999;
379  double WireLow = 999999;
380  fInterHigh_side = -999999;
381  fInterLow_side = 999999;
382  double min_ortdist(999), max_ortdist(-999); // needed to calculate width
383  //loop over all hits. Create coarse and fine profiles
384  // of the charge weight to help refine the start/end
385  // points and have a first guess of direction
386 
387  for (auto const& hit : fHitVector) {
388 
389  //calculate intercepts along
390  double intercept = hit.t - inv_2d_slope * (double)(hit.w);
391  double side_intercept = hit.t - fRough2DSlope * (double)(hit.w);
392 
393  util::PxPoint OnlinePoint;
394  // get coordinates of point on axis.
395  gser.GetPointOnLine(fRough2DSlope, fRough2DIntercept, &hit, OnlinePoint);
396 
397  double ortdist = gser.Get2DDistance(&OnlinePoint, &hit);
398  if (ortdist < min_ortdist) min_ortdist = ortdist;
399  if (ortdist > max_ortdist) max_ortdist = ortdist;
400 
401  if (inv_2d_slope != -999999) //almost all cases
402  {
403  if (intercept > InterHigh) { InterHigh = intercept; }
404 
405  if (intercept < InterLow) { InterLow = intercept; }
406  }
407  else //slope is practically horizontal. Care only about wires.
408  {
409  if (hit.w > WireHigh) { WireHigh = hit.w; }
410  if (hit.w < WireLow) { WireLow = hit.w; }
411  }
412 
413  if (side_intercept > fInterHigh_side) { fInterHigh_side = side_intercept; }
414 
415  if (side_intercept < fInterLow_side) { fInterLow_side = side_intercept; }
416  }
417  // end of first HitIter loop, at this point we should
418  // have the extreme intercepts
419 
420  /////////////////////////////////////////////
421  // Second loop. Fill profiles.
422  /////////////////////////////////////////////
423 
424  // get projected points at the beginning and end of the axis.
425 
426  util::PxPoint HighOnlinePoint, LowOnlinePoint, BeginOnlinePoint, EndOnlinePoint;
427 
428  if (inv_2d_slope != -999999) //almost all cases
429  {
430  gser.GetPointOnLineWSlopes(fRough2DSlope, fRough2DIntercept, InterHigh, HighOnlinePoint);
431  gser.GetPointOnLineWSlopes(fRough2DSlope, fRough2DIntercept, InterLow, LowOnlinePoint);
432  }
433  else //need better treatment of horizontal events.
434  {
435  util::PxPoint ptemphigh(fPlane, WireHigh, (fInterHigh_side + fInterLow_side) / 2);
436  util::PxPoint ptemplow(fPlane, WireLow, (fInterHigh_side + fInterLow_side) / 2);
437  gser.GetPointOnLine(fRough2DSlope, fRough2DIntercept, &ptemphigh, HighOnlinePoint);
438  gser.GetPointOnLine(fRough2DSlope, fRough2DIntercept, &ptemplow, LowOnlinePoint);
439  }
440  if (verbose) {
441  std::cout << " Low online point: " << LowOnlinePoint.w << ", " << LowOnlinePoint.t
442  << " High: " << HighOnlinePoint.w << ", " << HighOnlinePoint.t << std::endl;
443  //define BeginOnlinePoint as the one with lower wire number (for now), adjust intercepts accordingly
444  }
445  if (HighOnlinePoint.w >= LowOnlinePoint.w) {
446  BeginOnlinePoint = LowOnlinePoint;
447  fBeginIntercept = InterLow;
448  EndOnlinePoint = HighOnlinePoint;
449  fEndIntercept = InterHigh;
450  }
451  else {
452  BeginOnlinePoint = HighOnlinePoint;
453  fBeginIntercept = InterHigh;
454  EndOnlinePoint = LowOnlinePoint;
455  fEndIntercept = InterLow;
456  }
457 
458  fProjectedLength = gser.Get2DDistance(&BeginOnlinePoint, &EndOnlinePoint);
459 
460  /////////////////// THe binning is now set here
461  fCoarseNbins = 2;
462 
464  if (fProfileNbins < 10) fProfileNbins = 10;
465 
466  fChargeProfile.clear();
467  fCoarseChargeProfile.clear();
468  fChargeProfile.resize(fProfileNbins, 0);
470 
471  ////////////////////////// end of new binning
472  // Some fitting variables to make a histogram:
473 
474  // TODO this is nonsense for small clusters
475  std::vector<double> ort_profile;
476  const int NBINS = 100;
477  ort_profile.resize(NBINS);
478 
479  std::vector<double> ort_dist_vect;
480  ort_dist_vect.reserve(fHitVector.size());
481 
482  double current_maximum = 0;
483  for (auto& hit : fHitVector) {
484 
485  util::PxPoint OnlinePoint;
486  // get coordinates of point on axis.
487  gser.GetPointOnLine(fRough2DSlope, &BeginOnlinePoint, &hit, OnlinePoint);
488  double ortdist = gser.Get2DDistance(&OnlinePoint, &hit);
489 
490  double linedist = gser.Get2DDistance(&OnlinePoint, &BeginOnlinePoint);
491  ort_dist_vect.push_back(ortdist);
492  int ortbin;
493  if (ortdist == 0)
494  ortbin = 0;
495  else if (max_ortdist == min_ortdist)
496  ortbin = 0;
497  else
498  ortbin = (int)(ortdist - min_ortdist) / (max_ortdist - min_ortdist) * (NBINS - 1);
499 
500  ort_profile.at(ortbin) += hit.charge;
501 
502  //////////////////////////////////////////////////////////////////////
503  //calculate the weight along the axis, this formula is based on rough guessology.
504  // there is no physics motivation behind the particular numbers, A.S.
505  // \todo: switch to something that is motivated and easier to
506  // spell than guessology. C.A.
507  ///////////////////////////////////////////////////////////////////////
508  double weight = (ortdist < 1.) ? 10 * (hit.charge) : 5 * (hit.charge) / ortdist;
509 
510  int fine_bin =
511  fProjectedLength ? (int)(linedist / fProjectedLength * (fProfileNbins - 1)) : 0;
512  int coarse_bin =
513  fProjectedLength ? (int)(linedist / fProjectedLength * (fCoarseNbins - 1)) : 0;
514 
515  if (fine_bin < fProfileNbins) //only fill if bin number is in range
516  {
517  fChargeProfile.at(fine_bin) += weight;
518 
519  //find maximum bin on the fly:
520  if (fChargeProfile.at(fine_bin) > current_maximum && fine_bin != 0 &&
521  fine_bin != fProfileNbins - 1) {
522  current_maximum = fChargeProfile.at(fine_bin);
523  fProfileMaximumBin = fine_bin;
524  }
525  }
526 
527  if (coarse_bin < fCoarseNbins) //only fill if bin number is in range
528  fCoarseChargeProfile.at(coarse_bin) += weight;
529 
530  } // end second loop on hits. Now should have filled profile vectors.
531 
532  if (verbose) std::cout << "end second loop " << std::endl;
533 
534  double integral = 0;
535  int ix = 0;
536  for (ix = 0; ix < NBINS; ix++) {
537  integral += ort_profile.at(ix);
538  if (integral >= 0.95 * fParams.sum_charge) { break; }
539  }
540 
541  fParams.width =
542  2 * (double)ix / (double)(NBINS - 1) *
543  (double)(max_ortdist -
544  min_ortdist); // multiply by two because ortdist is folding in both sides.
545 
546  if (verbose) std::cout << " after width " << std::endl;
547 
548  if (drawOrtHistos) {
549  TH1F* ortDistHist = new TH1F(
550  "ortDistHist", "Orthogonal Distance to axis;Distance (cm)", 100, min_ortdist, max_ortdist);
551  TH1F* ortDistHistCharge =
552  new TH1F("ortDistHistCharge",
553  "Orthogonal Distance to axis (Charge Weighted);Distance (cm)",
554  100,
555  min_ortdist,
556  max_ortdist);
557  TH1F* ortDistHistAndrzej = new TH1F("ortDistHistAndrzej",
558  "Orthogonal Distance Andrzej weighting",
559  100,
560  min_ortdist,
561  max_ortdist);
562 
563  for (int index = 0; index < fParams.N_Hits; index++) {
564  ortDistHist->Fill(ort_dist_vect.at(index));
565  ortDistHistCharge->Fill(ort_dist_vect.at(index), fHitVector.at(index).charge);
566  double weight = (ort_dist_vect.at(index) < 1.) ?
567  10 * (fHitVector.at(index).charge) :
568  (5 * (fHitVector.at(index).charge) / ort_dist_vect.at(index));
569  ortDistHistAndrzej->Fill(ort_dist_vect.at(index), weight);
570  }
571  ortDistHist->Scale(1.0 / ortDistHist->Integral());
572  ortDistHistCharge->Scale(1.0 / ortDistHistCharge->Integral());
573  ortDistHistAndrzej->Scale(1.0 / ortDistHistAndrzej->Integral());
574 
575  TCanvas* ortCanv = new TCanvas("ortCanv", "ortCanv", 600, 600);
576  ortCanv->cd();
577  ortDistHistAndrzej->SetLineColor(3);
578  ortDistHistAndrzej->Draw("");
579  ortDistHistCharge->Draw("same");
580  ortDistHist->SetLineColor(2);
581  ortDistHist->Draw("same");
582 
583  TLegend* leg = new TLegend(.4, .6, .8, .9);
584  leg->SetHeader("Charge distribution");
585  leg->AddEntry(ortDistHist, "Unweighted Hits");
586  leg->AddEntry(ortDistHistCharge, "Charge Weighted Hits");
587  leg->AddEntry(ortDistHistAndrzej, "Charge Weighted by Guessology");
588  leg->Draw();
589 
590  ortCanv->Update();
591  }
592 
595 
596  //calculate the forward and backward integrals counting int the maximum bin.
597 
598  for (int ibin = 0; ibin < fProfileNbins; ibin++) {
601  }
602 
603  // now, we have the forward and backward integral so start
604  // stepping forward and backward to find the trunk of the
605  // shower. This is done making sure that the running
606  // integral until less than 1% is left and the bin is above
607  // a set theshold many empty bins.
608 
609  //forward loop
610  double running_integral = fProfileIntegralForward;
611  int startbin, endbin;
612  for (startbin = fProfileMaximumBin; startbin > 1 && startbin < (int)(fChargeProfile.size());
613  startbin--) {
614  running_integral -= fChargeProfile.at(startbin);
615  if (fChargeProfile.at(startbin) < fChargeCutoffThreshold.at(fPlane) &&
616  running_integral / fProfileIntegralForward < 0.01)
617  break;
618  else if (fChargeProfile.at(startbin) < fChargeCutoffThreshold.at(fPlane) &&
619  startbin - 1 > 0 &&
620  fChargeProfile.at(startbin - 1) < fChargeCutoffThreshold.at(fPlane) &&
621  startbin - 2 > 0 &&
622  fChargeProfile.at(startbin - 2) < fChargeCutoffThreshold.at(fPlane))
623  break;
624  }
625 
626  //backward loop
627  running_integral = fProfileIntegralBackward;
628  for (endbin = fProfileMaximumBin; endbin > 0 && endbin < fProfileNbins - 1; endbin++) {
629  running_integral -= fChargeProfile.at(endbin);
630  if (fChargeProfile.at(endbin) < fChargeCutoffThreshold.at(fPlane) &&
631  running_integral / fProfileIntegralBackward < 0.01)
632  break;
633  else if (fChargeProfile.at(endbin) < fChargeCutoffThreshold.at(fPlane) &&
634  endbin + 1 < fProfileNbins - 1 && endbin + 2 < fProfileNbins - 1 &&
635  fChargeProfile.at(endbin + 1) < fChargeCutoffThreshold.at(fPlane) &&
636  fChargeProfile.at(endbin + 2) < fChargeCutoffThreshold.at(fPlane))
637  break;
638  }
639 
640  // now have profile start and endpoints. Now translate to wire/time.
641  // Will use wire/time that are on the rough axis.
642  // fProjectedLength is the length from fInterLow to interhigh a
643  // long the rough_2d_axis on bin distance is:
644  // util::PxPoint OnlinePoint;
645 
646  if (inv_2d_slope != -999999) //almost all cases
647  {
648  double ort_intercept_begin =
649  fBeginIntercept + (fEndIntercept - fBeginIntercept) / fProfileNbins * startbin;
650  gser.GetPointOnLineWSlopes(
651  fRough2DSlope, fRough2DIntercept, ort_intercept_begin, fRoughBeginPoint);
652 
653  double ort_intercept_end =
654  fBeginIntercept + (fEndIntercept - fBeginIntercept) / fProfileNbins * endbin;
656 
657  gser.GetPointOnLineWSlopes(
658  fRough2DSlope, fRough2DIntercept, ort_intercept_end, fRoughEndPoint);
660  }
661  else {
662  double wire_begin = WireLow + (WireHigh - WireLow) / fProfileNbins * startbin;
663 
664  util::PxPoint ptemphigh(fPlane, wire_begin, (fInterHigh_side + fInterLow_side) / 2);
665  gser.GetPointOnLine(fRough2DSlope, fRough2DIntercept, &ptemphigh, fRoughBeginPoint);
667 
668  double wire_end = WireLow + (WireHigh - WireLow) / fProfileNbins * endbin;
669 
670  util::PxPoint ptemplow(fPlane, wire_end, (fInterHigh_side + fInterLow_side) / 2);
671  gser.GetPointOnLine(fRough2DSlope, fRough2DIntercept, &ptemplow, fRoughEndPoint);
673  }
674 
675  if (verbose)
676  std::cout << " rough start points " << fRoughBeginPoint.w << ", " << fRoughBeginPoint.t
677  << " end: " << fRoughEndPoint.w << " " << fRoughEndPoint.t << std::endl;
678 
679  // ok. now have fRoughBeginPoint and fRoughEndPoint. No decision about direction has been made yet.
682 
684 
685  fTimeRecord_ProcName.push_back("GetProfileInfo");
686  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
687  }
void GetRoughAxis(bool override=false)
std::vector< double > fCoarseChargeProfile
std::vector< std::string > fTimeRecord_ProcName
std::vector< util::PxHit > fHitVector
weight
Definition: test.py:257
double sum_charge
Sum charge of hits in ADC.
Definition: ClusterParams.h:27
cluster::cluster_params fParams
double t
Definition: PxUtils.h:11
std::vector< double > fTimeRecord_ProcTime
util::PxPoint start_point
start point
Definition: ClusterParams.h:24
Detector simulation of raw signals on wires.
double w
Definition: PxUtils.h:10
std::vector< double > fChargeProfile
util::PxPoint end_point
end point
Definition: ClusterParams.h:25
unsigned int plane
Definition: PxUtils.h:12
QTextStream & endl(QTextStream &s)
std::vector< double > fChargeCutoffThreshold
void cluster::ClusterParamsAlg::GetRoughAxis ( bool  override = false)

Calculates the following variables: verticalness fRough2DSlope fRough2DIntercept

Parameters
override[description]

Definition at line 279 of file ClusterParamsAlg.cxx.

280  {
281  if (!override) { //Override being set, we skip all this logic.
282  //OK, no override. Stop if we're already finshed.
283  if (fFinishedGetRoughAxis) return;
284  //Try to run the previous function if not yet done.
285  if (!fFinishedGetAverages) GetAverages(true);
286  }
287  else {
288  //Try to run the previous function if not yet done.
289  if (!fFinishedGetAverages) GetAverages(true);
290  }
291 
292  TStopwatch localWatch;
293  localWatch.Start();
294 
295  double rmsx = 0.0;
296  double rmsy = 0.0;
297  double rmsq = 0.0;
298  //using the charge weighted coordinates find the axis from slope
299  double ncw = 0;
300  double sumtime = 0; //from sum averages
301  double sumwire = 0; //from sum averages
302  double sumwiretime = 0; //sum over (wire*time)
303  double sumwirewire = 0; //sum over (wire*wire)
304  //next loop over all hits again
305 
306  for (auto& hit : fHitVector) {
307  // First, abuse this loop to calculate rms in x and y
308  rmsx += pow(fParams.mean_x - hit.w, 2) / fParams.N_Hits;
309  rmsy += pow(fParams.mean_y - hit.t, 2) / fParams.N_Hits;
310  rmsq += pow(fParams.mean_charge - hit.charge, 2) / fParams.N_Hits;
311  //if charge is above avg_charge
312 
313  if (hit.charge > fParams.mean_charge) {
314  ncw += 1;
315  sumwire += hit.w;
316  sumtime += hit.t;
317  sumwiretime += hit.w * hit.t;
318  sumwirewire += pow(hit.w, 2);
319  } //for high charge
320  } //For hh loop
321 
322  fParams.rms_x = sqrt(rmsx);
323  fParams.rms_y = sqrt(rmsy);
324  fParams.RMS_charge = sqrt(rmsq);
325 
326  fParams.N_Hits_HC = ncw;
327  //Looking for the slope and intercept of the line above avg_charge hits
328 
329  if ((ncw * sumwirewire - sumwire * sumwire) > 0.00001)
330  fRough2DSlope = (ncw * sumwiretime - sumwire * sumtime) /
331  (ncw * sumwirewire - sumwire * sumwire); //slope for cw
332  else
333  fRough2DSlope = 999;
336  fParams.charge_wgt_y - fRough2DSlope * (fParams.charge_wgt_x) : //intercept for cw
337  0.;
338 
339  //Getthe 2D_angle
340  fParams.cluster_angle_2d = atan(fRough2DSlope) * 180 / PI;
341 
342  fFinishedGetRoughAxis = true;
343 
344  fTimeRecord_ProcName.push_back("GetRoughAxis");
345  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
346  }
std::vector< std::string > fTimeRecord_ProcName
double charge_wgt_y
Mean of hits along y, charge weighted.
Definition: ClusterParams.h:38
double mean_x
Mean of hits along x, peaks only.
Definition: ClusterParams.h:33
constexpr T pow(T x)
Definition: pow.h:72
std::vector< util::PxHit > fHitVector
double mean_y
Mean of hits along y, peaks only.
Definition: ClusterParams.h:34
constexpr double PI
bool isfinite(Vector const &v)
Returns whether all components of the vector are finite.
cluster::cluster_params fParams
std::vector< double > fTimeRecord_ProcTime
void GetAverages(bool override=false)
Detector simulation of raw signals on wires.
double mean_charge
Mean (average) charge of hits in ADC.
Definition: ClusterParams.h:28
double charge_wgt_x
Mean of hits along x, charge weighted.
Definition: ClusterParams.h:37
double cluster_angle_2d
Linear best fit to high-charge hits in the cluster.
Definition: ClusterParams.h:39
double rms_x
rms of hits along x (wires)
Definition: ClusterParams.h:35
double rms_y
rms of hits along y, (time)
Definition: ClusterParams.h:36
void cluster::ClusterParamsAlg::Initialize ( void  )

Definition at line 133 of file ClusterParamsAlg.cxx.

134  {
135  fTimeRecord_ProcName.clear();
136  fTimeRecord_ProcTime.clear();
137 
138  TStopwatch localWatch;
139  localWatch.Start();
140 
141  //--- Initilize attributes values ---//
142  fFinishedGetAverages = false;
143  fFinishedGetRoughAxis = false;
144  fFinishedGetProfileInfo = false;
146  fFinishedRefineDirection = false;
147  fFinishedGetFinalSlope = false;
149  fFinishedTrackShowerSep = false;
150  fFinishedGetEndCharges = false;
151 
152  fRough2DSlope = -999.999; // slope
153  fRough2DIntercept = -999.999; // slope
154 
155  fRoughBeginPoint.w = -999.999;
156  fRoughEndPoint.w = -999.999;
157 
158  fRoughBeginPoint.t = -999.999;
159  fRoughEndPoint.t = -999.999;
160 
161  fProfileIntegralForward = 999.999;
162  fProfileIntegralBackward = 999.999;
163  fProfileMaximumBin = -999;
164 
165  fChargeCutoffThreshold.clear();
166  fChargeCutoffThreshold.reserve(3);
167  fChargeCutoffThreshold.push_back(500);
168  fChargeCutoffThreshold.push_back(500);
169  fChargeCutoffThreshold.push_back(1000);
170 
171  fHitVector.clear();
172 
173  fParams.Clear();
174 
175  fTimeRecord_ProcName.push_back("Initialize");
176  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
177  }
std::vector< std::string > fTimeRecord_ProcName
std::vector< util::PxHit > fHitVector
cluster::cluster_params fParams
double t
Definition: PxUtils.h:11
std::vector< double > fTimeRecord_ProcTime
double w
Definition: PxUtils.h:10
std::vector< double > fChargeCutoffThreshold
double cluster::ClusterParamsAlg::IntegrateFitCharge ( util::GeometryUtilities const &  gser,
double  from_length,
double  to_length,
unsigned int  fit_first_bin,
unsigned int  fit_end_bin 
)
protected

Integrates the charge between two positions in the cluster axis.

Parameters
from_lengthposition on the axis to start integration from, in cm
to_lengthposition on the axis to end integration at, in cm
fit_first_binfirst bin for the charge fit
fit_end_binnext-to-last bin for the charge fit
Returns
the charged fit integrated in the specified range, in ADC counts
See also
StartCharge(), EndCharge()

This function provides an almost-punctual charge at a position in the axis. Since the effective punctual charge is 0 ADC counts by definition, the charge can be integrated for some length. The procedure is made of two steps:

  1. the charge profile is parametrized with a linear fit within the specified region
  2. an integration of that fit is performed along the segment specified.

The region at point 1. is from fit_first_bin to fit_end_bin. These are specified in bin units. The binning is the one of the charge profile. It is suggested that a few bins are always kept, say 5 to 10, to reduce statistical fluctuations but maintaining a decent hypothesis of linearity along the range. The linear fit weighs all the bins in the profile the same.

The region at point to is from from_length to to_length, and it is measured in cm along the cluster axis, starting at the start of the cluster.

Definition at line 1408 of file ClusterParamsAlg.cxx.

1413  {
1414  // first compute the information on the charge profile
1415  GetProfileInfo(gser);
1416 
1417  // first things first
1418  if (fit_first_bin > fit_end_bin) std::swap(fit_first_bin, fit_end_bin);
1419 
1420  // how many bins can we use?
1421  const unsigned int nbins =
1422  std::min(fit_end_bin - fit_first_bin, (unsigned int)fChargeProfile.size());
1423  if (nbins == 0) return 0;
1424 
1425  // move the specified range within reason
1426  if (fit_end_bin > fChargeProfile.size()) {
1427  fit_end_bin = fChargeProfile.size();
1428  fit_first_bin = fit_end_bin - nbins;
1429  }
1430 
1431  // if we have to use only one bin, so be it
1432  if (nbins < 2) return fChargeProfile[fit_first_bin];
1433 
1434  // fit the charge profile vs. bin number;
1435  // we assume that the first bin (#0) starts from the very beginning of the
1436  // cluster, that is at axis coordinate 0.
1438  for (unsigned int iBin = fit_first_bin; iBin < fit_end_bin; ++iBin) {
1439  // should we use a Poisson error instead of no error?
1440  fitter.add((double)iBin, fChargeProfile[iBin]);
1441  } // for
1442 
1443  // get the linear fit parameters; [0] intercept [1] slope
1445  try {
1446  fit = fitter.FitParameters();
1447  }
1448  catch (std::range_error const&) { // oops...
1449  // this is actually unexpected, since the only reason for the fit to fail
1450  // would be a singular fit matrix, that should be pretty much impossible
1451  // given that the bin coordinates are well behaved and there are at least
1452  // two of them
1453  std::cerr << "IntegrateFitCharge(): linear fit failed!" << std::endl;
1454  return 0.;
1455  }
1456 
1457  // now determine the bins corresponding to the length to integrate;
1458  // note that length can be pathologic (0, negative...); not our problem!
1459  const double length_to_bin = (double(fChargeProfile.size() - 1) / fProjectedLength);
1460  const double from_bin = from_length * length_to_bin, to_bin = to_length * length_to_bin;
1461 
1462  // we want the integral between from_bin and to_bin now:
1463  return LinearIntegral(fit[1] /* m */, fit[0] /* q */, from_bin, to_bin);
1464  } // ClusterParamsAlg::IntegrateFitCharge()
void GetProfileInfo(util::GeometryUtilities const &gser, bool override=false)
Performs a linear regression of data.
Definition: SimpleFits.h:849
void swap(Handle< T > &a, Handle< T > &b)
std::array< Data_t, NParams > FitParameters_t
type of set of fit parameters
Definition: SimpleFits.h:866
std::vector< double > fChargeProfile
virtual FitParameters_t FitParameters() const override
Computes and returns all the parameters of the fit result.
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
bool add(Data_t x, Data_t y, Data_t sy=Data_t(1.0))
Definition: SimpleFits.h:366
static double LinearIntegral(double m, double q, double x1, double x2)
Returns the integral of f(x) = mx + q defined in [x1, x2].
QTextStream & endl(QTextStream &s)
double cluster::ClusterParamsAlg::LinearIntegral ( double  m,
double  q,
double  x1,
double  x2 
)
staticprotected

Returns the integral of f(x) = mx + q defined in [x1, x2].

Definition at line 1401 of file ClusterParamsAlg.cxx.

1402  {
1403  return m / 2. * (cet::square(x2) - cet::square(x1)) + q * (x2 - x1);
1404  }
constexpr T square(T x)
Definition: pow.h:21
size_t cluster::ClusterParamsAlg::MinNHits ( ) const
inline

Definition at line 38 of file ClusterParamsAlg.h.

39  {
40  return fMinNHits;
41  }
size_t fMinNHits
Cut value for # hits: below this value clusters are not evaluated.
float cluster::ClusterParamsAlg::MultipleHitDensity ( util::GeometryUtilities const &  gser)

Returns the number of multiple hits per wire.

Returns
the number of multiple hits per wire

This returns the number of wires with mmore than one hit belonging to this cluster, divided by the cluster length in cm.

Definition at line 1529 of file ClusterParamsAlg.cxx.

1530  {
1531  if (fHitVector.size() < 2) return 0.0F;
1532 
1533  // compute all the averages
1534  GetAverages();
1535  RefineStartPoints(gser); // fParams.length
1536 
1537  // return the relevant information
1538  return std::isnormal(fParams.length) ? fParams.multi_hit_wires / fParams.length : 0.;
1539  } // ClusterParamsAlg::MultipleHitDensity()
std::vector< util::PxHit > fHitVector
void RefineStartPoints(util::GeometryUtilities const &gser)
cluster::cluster_params fParams
void GetAverages(bool override=false)
float cluster::ClusterParamsAlg::MultipleHitWires ( )

Returns the number of multiple hits per wire.

Returns
the number of multiple hits per wire

This returns the fraction of wires that have more than one hit belonging to this cluster.

Definition at line 1516 of file ClusterParamsAlg.cxx.

1517  {
1518  if (fHitVector.size() < 2) return 0.0F;
1519 
1520  // compute all the averages
1521  GetAverages();
1522 
1523  // return the relevant information
1524  return std::isnormal(fParams.N_Wires) ? fParams.multi_hit_wires / fParams.N_Wires : 0.;
1525  } // ClusterParamsAlg::MultipleHitWires()
std::vector< util::PxHit > fHitVector
cluster::cluster_params fParams
void GetAverages(bool override=false)
int cluster::ClusterParamsAlg::Plane ( ) const
inline

Definition at line 284 of file ClusterParamsAlg.h.

285  {
286  return fPlane;
287  }
void cluster::ClusterParamsAlg::PrintFANNVector ( )

For debugging purposes, prints the result of GetFANNVector in a nicely formatted form.

Returns
[description]

Definition at line 110 of file ClusterParamsAlg.cxx.

111  {
112  std::vector<float> data;
113  GetFANNVector(data);
114  if (verbose) {
115  std::cout << "Printing FANN input vector:\n"
116  << " Opening Angle (normalized) ... : " << data[0] << "\n"
117  << " Opening Angle charge weight .. : " << data[1] << "\n"
118  << " Closing Angle (normalized) ... : " << data[2] << "\n"
119  << " Closing Angle charge weight .. : " << data[3] << "\n"
120  << " Principal Eigenvalue ......... : " << data[4] << "\n"
121  << " Secondary Eigenvalue ......... : " << data[5] << "\n"
122  << " Width / Length ............... : " << data[6] << "\n"
123  << " Hit Density / M.H.D. ......... : " << data[7] << "\n"
124  << " Percent MultiHit Wires ....... : " << data[8] << "\n"
125  << " Percent High Charge Hits ..... : " << data[9] << "\n"
126  << " Modified Hit Density ......... : " << data[10] << "\n"
127  << " Charge RMS / Mean Charge ...... : " << data[11] << "\n"
128  << " log(Sum Charge / Length) ...... : " << data[12] << "\n";
129  }
130  }
void GetFANNVector(std::vector< float > &data)
void cluster::ClusterParamsAlg::RefineDirection ( bool  override = false)

This function calculates the opening/closing angles It also makes a decision on whether or not to flip the direction the cluster. Then the start point is refined later.

Parameters
override[description]

Definition at line 1052 of file ClusterParamsAlg.cxx.

1053  {
1054  //
1055  // We don't use "override"? Should we remove? 05/01/14
1056  //
1057  if (!override) override = true;
1058 
1059  TStopwatch localWatch;
1060  localWatch.Start();
1061 
1062  // if(!override) { //Override being set, we skip all this logic.
1063  // //OK, no override. Stop if we're already finshed.
1064  // if (fFinishedRefineDirection) return;
1065  // //Try to run the previous function if not yet done.
1066  // if (!fFinishedGetProfileInfo) GetProfileInfo(true);
1067  // } else {
1068  // //Try to run the previous function if not yet done.
1069  // if (!fFinishedGetProfileInfo) GetProfileInfo(true);
1070  // }
1071 
1072  // double wire_2_cm = gser.WireToCm();
1073  // double time_2_cm = gser.TimeToCm();
1074 
1075  // Removing the conversion since these points are set above in cm-cm space
1076  // -- Corey
1077 
1078  util::PxPoint this_startPoint, this_endPoint;
1079 
1081  this_startPoint = fParams.start_point;
1082  this_endPoint = fParams.end_point;
1083  }
1084  else {
1085  this_startPoint = fRoughBeginPoint;
1086  this_endPoint = fRoughEndPoint;
1087  }
1088  if (verbose) {
1089  std::cout << "Angle: Start point: (" << this_startPoint.w << ", " << this_startPoint.t
1090  << ")\n";
1091  std::cout << "Angle: End point : (" << this_endPoint.w << ", " << this_endPoint.t << ")\n";
1092  }
1093  double endStartDiff_x = (this_endPoint.w - this_startPoint.w);
1094  double endStartDiff_y = (this_endPoint.t - this_startPoint.t);
1095  double rms_forward = 0;
1096  double rms_backward = 0;
1097  double mean_forward = 0;
1098  double mean_backward = 0;
1099  //double weight_total = 0;
1100  double hit_counter_forward = 0;
1101  double hit_counter_backward = 0;
1102 
1103  if (verbose && endStartDiff_y == 0 && endStartDiff_x == 0) {
1104  std::cerr << "Error: end point and start point are the same!\n";
1105 
1106  fTimeRecord_ProcName.push_back("RefineDirection");
1107  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1108  return;
1109  }
1110 
1111  double percentage = 0.90;
1112  double percentage_HC = 0.90 * fParams.N_Hits_HC / fParams.N_Hits;
1113  const int NBINS = 200;
1114  const double wgt = 1.0 / fParams.N_Hits;
1115 
1116  // Containers for the angle histograms
1117  std::vector<float> opening_angle_bin(NBINS, 0.0);
1118  std::vector<float> closing_angle_bin(NBINS, 0.0);
1119  std::vector<float> opening_angle_highcharge_bin(NBINS, 0.0);
1120  std::vector<float> closing_angle_highcharge_bin(NBINS, 0.0);
1121  std::vector<float> opening_angle_chargeWgt_bin(NBINS, 0.0);
1122  std::vector<float> closing_angle_chargeWgt_bin(NBINS, 0.0);
1123  //hard coding this for now, should use SetRefineDirectionQMin function
1124  fQMinRefDir = 25;
1125 
1126  int index = -1;
1127  for (auto& hit : fHitVector) {
1128  index++;
1129  //skip this hit if below minimum cutoff param
1130  if (hit.charge < fQMinRefDir) continue;
1131 
1132  //weight_total = hit.charge;
1133 
1134  // Compute forward mean
1135  double hitStartDiff_x = (hit.w - this_startPoint.w);
1136  double hitStartDiff_y = (hit.t - this_startPoint.t);
1137 
1138  if (hitStartDiff_x == 0 && hitStartDiff_y == 0) continue;
1139 
1140  double cosangle_start = (endStartDiff_x * hitStartDiff_x + endStartDiff_y * hitStartDiff_y);
1141  cosangle_start /= (pow(pow(endStartDiff_x, 2) + pow(endStartDiff_y, 2), 0.5) *
1142  pow(pow(hitStartDiff_x, 2) + pow(hitStartDiff_y, 2), 0.5));
1143 
1144  if (cosangle_start > 0) {
1145  // Only take into account for hits that are in "front"
1146  //no weighted average, works better as flat average w/ min charge cut
1147  mean_forward += cosangle_start;
1148  rms_forward += pow(cosangle_start, 2);
1149  hit_counter_forward++;
1150  }
1151 
1152  // Compute backward mean
1153  double hitEndDiff_x = (hit.w - this_endPoint.w);
1154  double hitEndDiff_y = (hit.t - this_endPoint.t);
1155  if (hitEndDiff_x == 0 && hitEndDiff_y == 0) continue;
1156 
1157  double cosangle_end = (endStartDiff_x * hitEndDiff_x + endStartDiff_y * hitEndDiff_y) * (-1.);
1158  cosangle_end /= (pow(pow(endStartDiff_x, 2) + pow(endStartDiff_y, 2), 0.5) *
1159  pow(pow(hitEndDiff_x, 2) + pow(hitEndDiff_y, 2), 0.5));
1160 
1161  if (cosangle_end > 0) {
1162  //no weighted average, works better as flat average w/ min charge cut
1163  mean_backward += cosangle_end;
1164  rms_backward += pow(cosangle_end, 2);
1165  hit_counter_backward++;
1166  }
1167 
1168  int N_bins_OPEN = (NBINS - 1) * acos(cosangle_start) / PI;
1169  int N_bins_CLOSE = (NBINS - 1) * acos(cosangle_end) / PI;
1170  if (N_bins_OPEN < 0) N_bins_OPEN = 0;
1171  if (N_bins_CLOSE < 0) N_bins_CLOSE = 0;
1172 
1173  opening_angle_chargeWgt_bin[N_bins_OPEN] += hit.charge / fParams.sum_charge;
1174  closing_angle_chargeWgt_bin[N_bins_CLOSE] += hit.charge / fParams.sum_charge;
1175  opening_angle_bin[N_bins_OPEN] += wgt;
1176  closing_angle_bin[N_bins_CLOSE] += wgt;
1177 
1178  //Also fill bins for high charge hits
1179  if (hit.charge > fParams.mean_charge) {
1180  opening_angle_highcharge_bin[N_bins_OPEN] += wgt;
1181  closing_angle_highcharge_bin[N_bins_CLOSE] += wgt;
1182  }
1183  }
1184 
1185  //initialize iterators for the angles
1186  int iBin(0), jBin(0), kBin(0), lBin(0), mBin(0), nBin(0);
1187  double percentage_OPEN(0.0), percentage_CLOSE(0.0), percentage_OPEN_HC(0.0),
1188  percentage_CLOSE_HC(0.0), //The last 2 are for High Charge (HC)
1189  percentage_OPEN_CHARGEWGT(0.0), percentage_CLOSE_CHARGEWGT(0.0);
1190 
1191  for (iBin = 0; percentage_OPEN <= percentage && iBin < NBINS; iBin++) {
1192  percentage_OPEN += opening_angle_bin[iBin];
1193  }
1194 
1195  for (jBin = 0; percentage_CLOSE <= percentage && jBin < NBINS; jBin++) {
1196  percentage_CLOSE += closing_angle_bin[jBin];
1197  }
1198 
1199  for (kBin = 0; percentage_OPEN_HC <= percentage_HC && kBin < NBINS; kBin++) {
1200  percentage_OPEN_HC += opening_angle_highcharge_bin[kBin];
1201  }
1202 
1203  for (lBin = 0; percentage_CLOSE_HC <= percentage_HC && lBin < NBINS; lBin++) {
1204  percentage_CLOSE_HC += closing_angle_highcharge_bin[lBin];
1205  }
1206 
1207  for (mBin = 0; percentage_OPEN_CHARGEWGT <= percentage && mBin < NBINS; mBin++) {
1208  percentage_OPEN_CHARGEWGT += opening_angle_chargeWgt_bin[mBin];
1209  }
1210 
1211  for (nBin = 0; percentage_CLOSE_CHARGEWGT <= percentage && nBin < NBINS; nBin++) {
1212  percentage_CLOSE_CHARGEWGT += closing_angle_chargeWgt_bin[nBin];
1213  }
1214 
1215  double opening_angle = iBin * PI / NBINS;
1216  double closing_angle = jBin * PI / NBINS;
1217  double opening_angle_highcharge = kBin * PI / NBINS;
1218  double closing_angle_highcharge = lBin * PI / NBINS;
1219  double opening_angle_charge_wgt = mBin * PI / NBINS;
1220  double closing_angle_charge_wgt = nBin * PI / NBINS;
1221 
1222  double value_1 = closing_angle / opening_angle - 1;
1223  double value_2 = (fCoarseChargeProfile[0] / fCoarseChargeProfile[1]);
1224  if (value_2 < 100.0 && value_2 > 0.01)
1225  value_2 = log(value_2);
1226  else
1227  value_2 = 0.0;
1228  double value_3 = closing_angle_charge_wgt / opening_angle_charge_wgt - 1;
1229 
1230  // Using a sigmoid function to determine flipping.
1231  // I am going to weigh all of the values above (1, 2, 3) equally.
1232  // But, since value 2 in particular, I'm going to cut things off at 5.
1233 
1234  if (value_1 > 3) value_1 = 3.0;
1235  if (value_1 < -3) value_1 = -3.0;
1236  if (value_2 > 3) value_2 = 3.0;
1237  if (value_2 < -3) value_2 = -3.0;
1238  if (value_3 > 3) value_3 = 3.0;
1239  if (value_3 < -3) value_3 = -3.0;
1240 
1241  double Exp = exp(value_1 + value_2 + value_3);
1242  double sigmoid = (Exp - 1) / (Exp + 1);
1243 
1244  bool flip = false;
1245  if (sigmoid < 0) flip = true;
1246  if (flip) {
1247  if (verbose) std::cout << "Flipping!" << std::endl;
1248  std::swap(opening_angle, closing_angle);
1249  std::swap(opening_angle_highcharge, closing_angle_highcharge);
1250  std::swap(opening_angle_charge_wgt, closing_angle_charge_wgt);
1253  }
1254  else if (verbose) {
1255  std::cout << "Not Flipping!\n";
1256  }
1257 
1258  fParams.opening_angle = opening_angle;
1259  fParams.opening_angle_charge_wgt = opening_angle_charge_wgt;
1260  fParams.closing_angle = closing_angle;
1261  fParams.closing_angle_charge_wgt = closing_angle_charge_wgt;
1262 
1263  fFinishedRefineDirection = true;
1264 
1265  fTimeRecord_ProcName.push_back("RefineDirection");
1266  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1267  } //end RefineDirection
double closing_angle_charge_wgt
Same for charge_wgt.
Definition: ClusterParams.h:44
std::vector< double > fCoarseChargeProfile
std::vector< std::string > fTimeRecord_ProcName
constexpr T pow(T x)
Definition: pow.h:72
std::vector< util::PxHit > fHitVector
constexpr double PI
double sum_charge
Sum charge of hits in ADC.
Definition: ClusterParams.h:27
cluster::cluster_params fParams
void swap(Handle< T > &a, Handle< T > &b)
double t
Definition: PxUtils.h:11
std::vector< double > fTimeRecord_ProcTime
util::PxPoint start_point
start point
Definition: ClusterParams.h:24
Detector simulation of raw signals on wires.
double w
Definition: PxUtils.h:10
util::PxPoint end_point
end point
Definition: ClusterParams.h:25
double mean_charge
Mean (average) charge of hits in ADC.
Definition: ClusterParams.h:28
double opening_angle_charge_wgt
Same for charge_wgt.
Definition: ClusterParams.h:42
double opening_angle
Width of angular distubtion wrt vertx.
Definition: ClusterParams.h:41
double closing_angle
Width of angular distubtion wrt endpoint.
Definition: ClusterParams.h:43
QTextStream & endl(QTextStream &s)
void cluster::ClusterParamsAlg::RefineStartPointAndDirection ( util::GeometryUtilities const &  gser,
bool  override = false 
)

Definition at line 1296 of file ClusterParamsAlg.cxx.

1297  {
1298  // This function is meant to pick the direction.
1299  // It refines both the start and end point, and then asks
1300  // if it should flip.
1301 
1302  TStopwatch localWatch;
1303  localWatch.Start();
1304 
1305  if (verbose) std::cout << " here!!! " << std::endl;
1306 
1307  if (!override) { //Override being set, we skip all this logic.
1308  //OK, no override. Stop if we're already finshed.
1310  fTimeRecord_ProcName.push_back("RefineStartPointAndDirection");
1311  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1312  return;
1313  }
1314  //Try to run the previous function if not yet done.
1315  if (!fFinishedGetProfileInfo) GetProfileInfo(gser, true);
1316  }
1317  else {
1318  //Try to run the previous function if not yet done.
1319  if (!fFinishedGetProfileInfo) GetProfileInfo(gser, true);
1320  }
1321  if (verbose) {
1322  std::cout << "REFINING .... " << std::endl;
1323  std::cout << " Rough start and end point: " << std::endl;
1324  std::cout << " s: (" << fParams.start_point.w << ", " << fParams.start_point.t << ")"
1325  << std::endl;
1326  std::cout << " e: (" << fParams.end_point.w << ", " << fParams.end_point.t << ")"
1327  << std::endl;
1328  }
1329  RefineStartPoints(gser);
1330  if (verbose) {
1331  std::cout << " Once Refined start and end point: " << std::endl;
1332  std::cout << " s: (" << fParams.start_point.w << ", " << fParams.start_point.t << ")"
1333  << std::endl;
1334  std::cout << " e: (" << fParams.end_point.w << ", " << fParams.end_point.t << ")"
1335  << std::endl;
1338  }
1339  RefineStartPoints(gser);
1340  if (verbose) {
1341  std::cout << " Twice Refined start and end point: " << std::endl;
1342  std::cout << " s: (" << fParams.start_point.w << ", " << fParams.start_point.t << ")"
1343  << std::endl;
1344  std::cout << " e: (" << fParams.end_point.w << ", " << fParams.end_point.t << ")"
1345  << std::endl;
1348  }
1349  RefineDirection();
1350  if (verbose) {
1351  std::cout << " Final start and end point: " << std::endl;
1352  std::cout << " s: (" << fParams.start_point.w << ", " << fParams.start_point.t << ")"
1353  << std::endl;
1354  std::cout << " e: (" << fParams.end_point.w << ", " << fParams.end_point.t << ")"
1355  << std::endl;
1356  }
1358 
1359  fTimeRecord_ProcName.push_back("RefineStartPointAndDirection");
1360  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1361  }
std::vector< std::string > fTimeRecord_ProcName
void GetProfileInfo(util::GeometryUtilities const &gser, bool override=false)
void RefineStartPoints(util::GeometryUtilities const &gser)
cluster::cluster_params fParams
void swap(Handle< T > &a, Handle< T > &b)
double t
Definition: PxUtils.h:11
std::vector< double > fTimeRecord_ProcTime
util::PxPoint start_point
start point
Definition: ClusterParams.h:24
double w
Definition: PxUtils.h:10
void RefineDirection(bool override=false)
util::PxPoint end_point
end point
Definition: ClusterParams.h:25
QTextStream & endl(QTextStream &s)
void cluster::ClusterParamsAlg::RefineStartPoints ( util::GeometryUtilities const &  gser)

Calculates the following variables: length width

Calculates the following variables: length width hit_density_1D hit_density_2D direction

Parameters
override[description]

Definition at line 699 of file ClusterParamsAlg.cxx.

700  {
701  TStopwatch localWatch;
702  localWatch.Start();
703 
704  // need to define physical direction with openind angles and pass that to Ryan's line finder.
705 
706  // Direction decision has been moved entirely to Refine Direction()
707  // opening and closing angles are already set
708  // they are associated with the start and endpoints.
709  // If refine direction opted to switch the shower direction,
710  // it also switched opening and closing angles.
711 
712  /////////////////////////////////
713  // fParams.direction= ...
714 
715  // Direction is decided by RefineDirection()
716  util::PxPoint startHit, endHit;
717  startHit = fRoughBeginPoint;
718  endHit = fRoughEndPoint;
719 
720  ///////////////////////////////
721  //Ryan's Shower Strip finder work here.
722  //First we need to define the strip width that we want
723  double d = 0.6; //this is the width of the strip.... this needs to be tuned to something.
724  //===============================================================================================================
725  // Will need to feed in the set of hits that we want.
726  // const std::vector<util::PxHit*> whole;
727  std::vector<const util::PxHit*> subhit;
728  double linearlimit = 8;
729  double ortlimit = 12;
730  double lineslopetest(fRough2DSlope);
731  util::PxHit averageHit;
732  //also are we sure this is actually doing what it is supposed to???
733  //are we sure this works?
734  //is anybody even listening? Were they eaten by a grue?
735  gser.SelectLocalHitlist(
736  fHitVector, subhit, startHit, linearlimit, ortlimit, lineslopetest, averageHit);
737 
738  if (!(subhit.size()) || subhit.size() < 3) {
739  if (verbose)
740  std::cout << "Subhit list is empty or too small. Using rough start/end points..."
741  << std::endl;
744 
746 
747  fTimeRecord_ProcName.push_back("RefineStartPoints");
748  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
749 
750  return;
751  }
752 
753  double avgwire = averageHit.w;
754  double avgtime = averageHit.t;
755  //vertex in tilda-space pair(x-til,y-til)
756  std::vector<std::pair<double, double>> vertil;
757  //vector of the sum of the distance of a vector to every vertex in tilda-space
758  std::vector<double> vs;
759  // $$This needs to be corrected//this is the good hits that are between strip
760  std::vector<const util::PxHit*> ghits;
761  ghits.reserve(subhit.size());
762  int n = 0;
763  double fardistcurrent = 0;
764  util::PxHit startpoint;
765  double gwiretime = 0;
766  double gwire = 0;
767  double gtime = 0;
768  double gwirewire = 0;
769  //KAZU!!! I NEED this too//this will need to come from somewhere...
770  //This is some hit that is from the way far side of the entire shower cluster...
771  util::PxPoint farhit = fRoughEndPoint;
772 
773  //=============building the local vector===================
774  //this is clumsy... but just want to get something to work for now.
775  //THIS is REAL Horrible and CLUMSY... We can make things into helper functions soon.
776  std::vector<util::PxHit> returnhits;
777  std::vector<double> radiusofhit;
778  std::vector<int> closehits;
779  //==============================================================================
780  //Now we need to do the transformation into "tilda-space"
781  for (unsigned int a = 0; a < subhit.size(); a++) {
782  for (unsigned int b = a + 1; b < subhit.size(); b++) {
783  if (subhit.at(a)->w != subhit.at(b)->w) {
784  double xtil = ((subhit.at(a)->t - avgtime) - (subhit.at(b)->t - avgtime));
785  xtil /= ((subhit.at(a)->w - avgwire) - (subhit.at(b)->w - avgwire));
786  double ytil = (subhit.at(a)->w - avgwire) * xtil - (subhit.at(a)->t - avgtime);
787  //push back the tilda vertex point on the pair
788  std::pair<double, double> tv(xtil, ytil);
789  vertil.push_back(tv);
790  } //if Wires are not the same
791  } //for over b
792  } //for over a
793  // look at the distance from a tilda-vertex to all other tilda-verticies
794  for (unsigned int z = 0; z < vertil.size(); z++) {
795  double vd = 0; //the distance for vertex... just needs to be something 0
796  for (unsigned int c = 0; c < vertil.size(); c++)
797  vd += std::hypot(vertil.at(z).first - vertil.at(c).first,
798  vertil.at(z).second - vertil.at(c).second);
799  vs.push_back(vd);
800  vd = 0.0;
801  } //for z loop
802 
803  if (vs.size() == 0) //al hits on same wire?!
804  {
805  if (verbose)
806  std::cout << "vertil list is empty. all subhits are on the same wire?" << std::endl;
809 
811 
812  fTimeRecord_ProcName.push_back("RefineStartPoints");
813  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
814  return;
815  }
816  //need to find the min of the distance of vertex in tilda-space
817  //this will get me the area where things are most linear
818  int minvs = std::min_element(vs.begin(), vs.end()) - vs.begin();
819  // now use the min position to find the vertex in tilda-space
820  //now need to look a which hits are between the tilda lines from the points
821  //in the tilda space everything in wire time is based on the new origin which is at the average wire/time
822  double tilwire = vertil.at(minvs).first; //slope
823  double tiltimet = -vertil.at(minvs).second +
824  d * sqrt(1 + pow(tilwire, 2)); //negative cept is accounted for top strip
825  double tiltimeb = -vertil.at(minvs).second -
826  d * sqrt(1 + pow(tilwire, 2)); //negative cept is accounted for bottom strip
827  // look over the subhit list and ask for which are inside of the strip
828  for (unsigned int a = 0; a < subhit.size(); a++) {
829  double dtstrip =
830  (-tilwire * (subhit.at(a)->w - avgwire) + (subhit.at(a)->t - avgtime) - tiltimet) /
831  sqrt(tilwire * tilwire + 1);
832  double dbstrip =
833  (-tilwire * (subhit.at(a)->w - avgwire) + (subhit.at(a)->t - avgtime) - tiltimeb) /
834  sqrt(tilwire * tilwire + 1);
835 
836  if ((dtstrip < 0.0 && dbstrip > 0.0) || (dbstrip < 0.0 && dtstrip > 0.0)) {
837  ghits.push_back(subhit.at(a));
838  } //if between strip
839  } //for a loop
840  //=======================Do stuff with the good hits============================
841 
842  //of these small set of hits just fit a simple line.
843  //then we will need to see which of these hits is farthest away
844 
845  for (unsigned int g = 0; g < ghits.size(); g++) {
846  // should call the helper funtion to do the fit
847  //but for now since there is no helper function I will just write it here......again
848  n += 1;
849  gwiretime += ghits.at(g)->w * ghits.at(g)->t;
850  gwire += ghits.at(g)->w;
851  gtime += ghits.at(g)->t;
852  gwirewire += ghits.at(g)->w * ghits.at(g)->w;
853  // now work on calculating the distance in wire time space from the far point
854  //farhit needs to be a hit that is given to me
855  double fardist = std::hypot(ghits.at(g)->w - farhit.w, ghits.at(g)->t - farhit.t);
856  //come back to this... there is a better way to do this probably in the loop
857  //there should also be a check that the hit that is farthest away has subsequent hits after it on a few wires
858  if (fardist > fardistcurrent) {
859  fardistcurrent = fardist;
860  //if fardist... this is the point to use for the start point
861  startpoint.t = ghits.at(g)->t;
862  startpoint.w = ghits.at(g)->w;
863  startpoint.plane = ghits.at(g)->plane;
864  startpoint.charge = ghits.at(g)->charge;
865  }
866  } //for ghits loop
867  //This can be the new start point
868 
869  //should this be here? Id argue this might be done once outside.
870  fParams.length = gser.Get2DDistance((util::PxPoint*)&startpoint, &endHit);
871  if (fParams.length)
873  else
875 
876  if (fParams.length && fParams.width)
878  else
880 
881  fParams.start_point = startpoint;
882 
884 
885  fTimeRecord_ProcName.push_back("RefineStartPoints");
886  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
887  }
std::vector< std::string > fTimeRecord_ProcName
static constexpr double g
Definition: Units.h:144
constexpr T pow(T x)
Definition: pow.h:72
std::vector< util::PxHit > fHitVector
std::enable_if_t< std::is_arithmetic_v< T >, T > hypot(T x, T y)
Definition: hypot.h:60
cluster::cluster_params fParams
std::void_t< T > n
const double a
double t
Definition: PxUtils.h:11
std::vector< double > fTimeRecord_ProcTime
util::PxPoint start_point
start point
Definition: ClusterParams.h:24
double w
Definition: PxUtils.h:10
util::PxPoint end_point
end point
Definition: ClusterParams.h:25
double charge
area charge
Definition: PxUtils.h:39
static bool * b
Definition: config.cpp:1043
unsigned int plane
Definition: PxUtils.h:12
QTextStream & endl(QTextStream &s)
const util::PxPoint& cluster::ClusterParamsAlg::RoughEndPoint ( )
inline

Definition at line 187 of file ClusterParamsAlg.h.

188  {
189  return fRoughEndPoint;
190  }
double cluster::ClusterParamsAlg::RoughIntercept ( )
inline

Definition at line 198 of file ClusterParamsAlg.h.

199  {
200  return fRough2DIntercept;
201  }
double cluster::ClusterParamsAlg::RoughSlope ( )
inline

Definition at line 193 of file ClusterParamsAlg.h.

194  {
195  return fRough2DSlope;
196  }
const util::PxPoint& cluster::ClusterParamsAlg::RoughStartPoint ( )
inline

Definition at line 182 of file ClusterParamsAlg.h.

183  {
184  return fRoughBeginPoint;
185  }
int cluster::ClusterParamsAlg::SetHits ( const std::vector< util::PxHit > &  inhitlist)

Definition at line 51 of file ClusterParamsAlg.cxx.

52  {
53  Initialize();
54 
55  // Make default values
56  // Is done by the struct
57  if (!(inhitlist.size())) {
58  throw CRUException("Provided empty hit list!");
59  return -1;
60  }
61 
62  fHitVector.reserve(inhitlist.size());
63  for (auto h : inhitlist)
64  fHitVector.push_back(h);
65 
66  fPlane = fHitVector[0].plane;
67 
68  if (fHitVector.size() < fMinNHits) {
69  if (verbose)
70  std::cout << " the hitlist is too small. Continuing to run may result in crash!!! "
71  << std::endl;
72  return -1;
73  }
74 
75  return fHitVector.size();
76  }
std::vector< util::PxHit > fHitVector
size_t fMinNHits
Cut value for # hits: below this value clusters are not evaluated.
QTextStream & endl(QTextStream &s)
void cluster::ClusterParamsAlg::SetMinNHits ( size_t  nhit)
inline

Definition at line 32 of file ClusterParamsAlg.h.

33  {
34  fMinNHits = nhit;
35  }
uint8_t nhit
Definition: CRTFragment.hh:201
size_t fMinNHits
Cut value for # hits: below this value clusters are not evaluated.
void cluster::ClusterParamsAlg::setNeuralNetPath ( std::string  s)
inline

Definition at line 172 of file ClusterParamsAlg.h.

173  {
174  fNeuralNetPath = s;
175  }
static QCString * s
Definition: config.cpp:1042
void cluster::ClusterParamsAlg::SetPlane ( int  p)

Definition at line 79 of file ClusterParamsAlg.cxx.

80  {
81  fPlane = p;
82  for (auto& h : fHitVector)
83  h.plane = p;
86  }
std::vector< util::PxHit > fHitVector
cluster::cluster_params fParams
p
Definition: test.py:223
util::PxPoint start_point
start point
Definition: ClusterParams.h:24
util::PxPoint end_point
end point
Definition: ClusterParams.h:25
unsigned int plane
Definition: PxUtils.h:12
void cluster::ClusterParamsAlg::SetRefineDirectionQMin ( double  qmin)
inline

Definition at line 46 of file ClusterParamsAlg.h.

47  {
48  fQMinRefDir = qmin;
49  }
void cluster::ClusterParamsAlg::SetVerbose ( bool  yes = true)
inline

Definition at line 52 of file ClusterParamsAlg.h.

double cluster::ClusterParamsAlg::StartCharge ( util::GeometryUtilities const &  gser,
float  length = 1.,
unsigned int  nbins = 10 
)

Returns the expected charge at the beginning of the cluster.

Parameters
nbinsuse at least this number of charge bins from charge profile
lengthspace at the start of cluster where to collect charge, in cm the expected charge at the beginning of the cluster
See also
EndCharge(), IntegrateFitCharge()

ClusterParamsAlg extracts a binned charge profile, parametrized versus the distance from the start of the cluster. All the charge on the plane orthogonal to cluster axis is collapsed into the point where that plane intersects the axis. The resulting 1D distribution is then binned.

This method returns the charge under the first length cm of the cluster.

This method considers the first nbins of this charge distribution and through a linear fit determines the expected charge at the first bin. Then, it scales the result to reflect how much charge would be deposited in a space of length centimetres, according to this linear fit.

Note that length may be 0 (charge will be 0) or negative (sort of extrapolation ahead of the cluster start).

For more details, see IntegrateFitCharge().

Definition at line 1468 of file ClusterParamsAlg.cxx.

1471  {
1472  switch (fHitVector.size()) {
1473  case 0: return 0.;
1474  case 1:
1475  return fHitVector.back().charge;
1476  // the "default" is the rest of the function
1477  } // switch
1478 
1479  // this is the range of the fit:
1480  const unsigned int fit_first_bin = 0, fit_last_bin = nbins;
1481 
1482  return IntegrateFitCharge(gser, 0., length, fit_first_bin, fit_last_bin);
1483  } // ClusterParamsAlg::StartCharge()
double IntegrateFitCharge(util::GeometryUtilities const &gser, double from_length, double to_length, unsigned int fit_first_bin, unsigned int fit_end_bin)
Integrates the charge between two positions in the cluster axis.
std::vector< util::PxHit > fHitVector
template<typename Stream >
void cluster::ClusterParamsAlg::TimeReport ( Stream &  stream) const

Definition at line 402 of file ClusterParamsAlg.h.

403  {
404 
405  stream << " <<ClusterParamsAlg::TimeReport>> starts...\n";
406  for (size_t i = 0; i < fTimeRecord_ProcName.size(); ++i) {
407 
408  stream << " Function: " << fTimeRecord_ProcName[i].c_str()
409  << " ... Time = " << fTimeRecord_ProcTime[i] << " [s]\n";
410  }
411  stream << " <<ClusterParamsAlg::TimeReport>> ends...\n";
412  }
std::vector< std::string > fTimeRecord_ProcName
std::vector< double > fTimeRecord_ProcTime
void cluster::ClusterParamsAlg::TrackShowerSeparation ( bool  override = false)

Definition at line 1364 of file ClusterParamsAlg.cxx.

1365  {
1366  if (!override) return;
1367  if (verbose) std::cout << " ---- Trying T/S sep. ------ \n";
1368  if (enableFANN) {
1369  if (verbose) std::cout << " ---- Doing T/S sep. ------- \n";
1370  std::vector<float> FeatureVector, outputVector;
1371  GetFANNVector(FeatureVector);
1372  }
1373  else {
1374  if (verbose) std::cout << " ---- Failed T/S sep. ------ \n";
1375  }
1376  }
void GetFANNVector(std::vector< float > &data)

Member Data Documentation

bool cluster::ClusterParamsAlg::enableFANN
protected

Definition at line 343 of file ClusterParamsAlg.h.

double cluster::ClusterParamsAlg::fBeginIntercept
protected

Definition at line 323 of file ClusterParamsAlg.h.

std::vector<double> cluster::ClusterParamsAlg::fChargeCutoffThreshold
protected

Definition at line 304 of file ClusterParamsAlg.h.

std::vector<double> cluster::ClusterParamsAlg::fChargeProfile
protected

Definition at line 310 of file ClusterParamsAlg.h.

std::vector<double> cluster::ClusterParamsAlg::fChargeProfileNew
protected

Definition at line 313 of file ClusterParamsAlg.h.

std::vector<double> cluster::ClusterParamsAlg::fCoarseChargeProfile
protected

Definition at line 311 of file ClusterParamsAlg.h.

int cluster::ClusterParamsAlg::fCoarseNbins
protected

Definition at line 315 of file ClusterParamsAlg.h.

double cluster::ClusterParamsAlg::fEndIntercept
protected

Definition at line 324 of file ClusterParamsAlg.h.

bool cluster::ClusterParamsAlg::fFinishedGetAverages
protected

Definition at line 329 of file ClusterParamsAlg.h.

bool cluster::ClusterParamsAlg::fFinishedGetEndCharges
protected

Definition at line 337 of file ClusterParamsAlg.h.

bool cluster::ClusterParamsAlg::fFinishedGetFinalSlope
protected

Definition at line 334 of file ClusterParamsAlg.h.

bool cluster::ClusterParamsAlg::fFinishedGetProfileInfo
protected

Definition at line 331 of file ClusterParamsAlg.h.

bool cluster::ClusterParamsAlg::fFinishedGetRoughAxis
protected

Definition at line 330 of file ClusterParamsAlg.h.

bool cluster::ClusterParamsAlg::fFinishedRefineDirection
protected

Definition at line 333 of file ClusterParamsAlg.h.

bool cluster::ClusterParamsAlg::fFinishedRefineStartPointAndDirection
protected

Definition at line 335 of file ClusterParamsAlg.h.

bool cluster::ClusterParamsAlg::fFinishedRefineStartPoints
protected

Definition at line 332 of file ClusterParamsAlg.h.

bool cluster::ClusterParamsAlg::fFinishedTrackShowerSep
protected

Definition at line 336 of file ClusterParamsAlg.h.

std::vector<util::PxHit> cluster::ClusterParamsAlg::fHitVector
protected

This vector holds the pointer to hits. This should be used for computation for speed.

Definition at line 298 of file ClusterParamsAlg.h.

double cluster::ClusterParamsAlg::fInterHigh_side
protected

Definition at line 325 of file ClusterParamsAlg.h.

double cluster::ClusterParamsAlg::fInterLow_side
protected

Definition at line 326 of file ClusterParamsAlg.h.

size_t cluster::ClusterParamsAlg::fMinNHits
protected

Cut value for # hits: below this value clusters are not evaluated.

Definition at line 292 of file ClusterParamsAlg.h.

std::string cluster::ClusterParamsAlg::fNeuralNetPath

Definition at line 385 of file ClusterParamsAlg.h.

cluster::cluster_params cluster::ClusterParamsAlg::fParams

Definition at line 383 of file ClusterParamsAlg.h.

int cluster::ClusterParamsAlg::fPlane
protected

Definition at line 305 of file ClusterParamsAlg.h.

double cluster::ClusterParamsAlg::fProfileIntegralBackward
protected

Definition at line 319 of file ClusterParamsAlg.h.

double cluster::ClusterParamsAlg::fProfileIntegralForward
protected

Definition at line 318 of file ClusterParamsAlg.h.

int cluster::ClusterParamsAlg::fProfileMaximumBin
protected

Definition at line 317 of file ClusterParamsAlg.h.

int cluster::ClusterParamsAlg::fProfileNbins
protected

Definition at line 316 of file ClusterParamsAlg.h.

double cluster::ClusterParamsAlg::fProjectedLength
protected

Definition at line 320 of file ClusterParamsAlg.h.

double cluster::ClusterParamsAlg::fQMinRefDir
protected

Definition at line 308 of file ClusterParamsAlg.h.

double cluster::ClusterParamsAlg::fRough2DIntercept
protected

Definition at line 340 of file ClusterParamsAlg.h.

double cluster::ClusterParamsAlg::fRough2DSlope
protected

Definition at line 339 of file ClusterParamsAlg.h.

util::PxPoint cluster::ClusterParamsAlg::fRoughBeginPoint
protected

Definition at line 341 of file ClusterParamsAlg.h.

util::PxPoint cluster::ClusterParamsAlg::fRoughEndPoint
protected

Definition at line 342 of file ClusterParamsAlg.h.

std::vector<std::string> cluster::ClusterParamsAlg::fTimeRecord_ProcName

Definition at line 387 of file ClusterParamsAlg.h.

std::vector<double> cluster::ClusterParamsAlg::fTimeRecord_ProcTime

Definition at line 388 of file ClusterParamsAlg.h.

bool cluster::ClusterParamsAlg::verbose
protected

Definition at line 301 of file ClusterParamsAlg.h.


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