LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
cluster::ClusterParamsAlg Class Reference

#include "ClusterParamsAlg.h"

Public Member Functions

 ClusterParamsAlg ()
 Default constructor. More...
 
 ClusterParamsAlg (const std::vector< util::PxHit > &)
 Alternative constructor with larutil::PxHit vector. More...
 
 ~ClusterParamsAlg ()
 
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 Report (Stream &stream) const
 
template<typename Stream >
void TimeReport (Stream &stream) const
 
void GetFANNVector (std::vector< float > &data)
 
void PrintFANNVector ()
 
void FillParams (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 (bool override=false)
 
void RefineStartPoints (bool override=false)
 
void GetFinalSlope (bool override=false)
 
void GetEndCharges (bool override_=false)
 
void RefineDirection (bool override=false)
 
void RefineStartPointAndDirection (bool override=false)
 
void TrackShowerSeparation (bool override=false)
 
void setNeuralNetPath (std::string s)
 
void FillPolygon ()
 
void GetOpeningAngle ()
 
const util::PxPointRoughStartPoint ()
 
const util::PxPointRoughEndPoint ()
 
double RoughSlope ()
 
double RoughIntercept ()
 
double StartCharge (float length=1., unsigned int nbins=10)
 Returns the expected charge at the beginning of the cluster. More...
 
double EndCharge (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 ()
 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::string > fTimeRecord_ProcName
 
std::vector< double > fTimeRecord_ProcTime
 

Protected Member Functions

double IntegrateFitCharge (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

util::GeometryUtilitiesfGSer
 
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 42 of file ClusterParamsAlg.h.

Constructor & Destructor Documentation

cluster::ClusterParamsAlg::ClusterParamsAlg ( )

Default constructor.

Definition at line 25 of file ClusterParamsAlg.cxx.

References enableFANN, fGSer, fMinNHits, Initialize(), and verbose.

26  {
27  fMinNHits = 10;
28  fGSer=nullptr;
29  enableFANN = false;
30  verbose=true;
31  Initialize();
32  }
util::GeometryUtilities * fGSer
size_t fMinNHits
Cut value for # hits: below this value clusters are not evaluated.
cluster::ClusterParamsAlg::ClusterParamsAlg ( const std::vector< util::PxHit > &  inhitlist)

Alternative constructor with larutil::PxHit vector.

Definition at line 43 of file ClusterParamsAlg.cxx.

References enableFANN, fGSer, fMinNHits, SetHits(), and verbose.

44  {
45  fMinNHits = 10;
46  fGSer=nullptr;
47  enableFANN = false;
48  verbose=true;
49  SetHits(inhitlist);
50  }
util::GeometryUtilities * fGSer
size_t fMinNHits
Cut value for # hits: below this value clusters are not evaluated.
int SetHits(const std::vector< util::PxHit > &)
cluster::ClusterParamsAlg::~ClusterParamsAlg ( )
inline

Definition at line 55 of file ClusterParamsAlg.h.

References Initialize().

55 {}

Member Function Documentation

void cluster::ClusterParamsAlg::DisableFANN ( )
inline

Definition at line 268 of file ClusterParamsAlg.h.

References enableFANN.

void cluster::ClusterParamsAlg::EnableFANN ( )

Definition at line 236 of file ClusterParamsAlg.cxx.

References enableFANN.

Referenced by RoughIntercept().

236  {
237  enableFANN = true;
238  // ::cluster::FANNService::GetME()->GetFANNModule().LoadFromFile(fNeuralNetPath);
239  return;
240  }
double cluster::ClusterParamsAlg::EndCharge ( 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 1787 of file ClusterParamsAlg.cxx.

References fChargeProfile, fHitVector, fProjectedLength, GetProfileInfo(), and IntegrateFitCharge().

Referenced by cluster::StandardClusterParamsAlg::EndCharge(), GetEndCharges(), RoughIntercept(), and StartCharge().

1788  {
1789  switch (fHitVector.size()) {
1790  case 0: return 0.;
1791  case 1: return fHitVector.back().charge;
1792  // the "default" is the rest of the function
1793  } // switch
1794 
1795 
1796  // need the available number of bins and the axis length
1797  GetProfileInfo();
1798  const unsigned int MaxBins = fChargeProfile.size();
1799 
1800  // this is the range of the fit:
1801  const unsigned int fit_first_bin = MaxBins > nbins? MaxBins - nbins: 0,
1802  fit_last_bin = MaxBins;
1803 
1804  // now determine the integration range, in bin units;
1805  // get to the end, and go length backward;
1806  // note that length can be pathologic (0, negative...); not our problem!
1807  const double from = fProjectedLength -length, to = fProjectedLength;
1808 
1809  return IntegrateFitCharge(from, to, fit_first_bin, fit_last_bin);
1810  } // ClusterParamsAlg::EndCharge()
void GetProfileInfo(bool override=false)
std::vector< util::PxHit > fHitVector
double IntegrateFitCharge(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< double > fChargeProfile
void cluster::ClusterParamsAlg::FillParams ( 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 258 of file ClusterParamsAlg.cxx.

References GetAverages(), GetEndCharges(), GetFinalSlope(), GetProfileInfo(), GetRoughAxis(), RefineStartPointAndDirection(), and TrackShowerSeparation().

Referenced by SetVerbose().

266  {
267  GetAverages (override_DoGetAverages );
268  GetRoughAxis (override_DoGetRoughAxis );
269  GetProfileInfo (override_DoGetProfileInfo );
270  RefineStartPointAndDirection(override_DoStartPointsAndDirection);
271  // RefineDirection (override_DoRefineDirection );
272  // RefineStartPoints(override_DoRefineStartPoints);
273  GetFinalSlope (override_DoGetFinalSlope );
274  GetEndCharges (override_DoEndCharge);
275  TrackShowerSeparation(override_DoTrackShowerSep);
276  }
void GetRoughAxis(bool override=false)
void GetProfileInfo(bool override=false)
void GetFinalSlope(bool override=false)
void TrackShowerSeparation(bool override=false)
void GetAverages(bool override=false)
void RefineStartPointAndDirection(bool override=false)
void GetEndCharges(bool override_=false)
void cluster::ClusterParamsAlg::FillPolygon ( )

Definition at line 1568 of file ClusterParamsAlg.cxx.

References fGSer, fHitVector, fParams, fTimeRecord_ProcName, fTimeRecord_ProcTime, cluster::cluster_params::PolyObject, and util::GeometryUtilities::SelectPolygonHitList().

Referenced by setNeuralNetPath().

1569  {
1570 
1571  TStopwatch localWatch;
1572  localWatch.Start();
1573 
1574  if(fHitVector.size()) {
1575  std::vector<const util::PxHit*> container_polygon;
1576  fGSer->SelectPolygonHitList(fHitVector,container_polygon);
1577  //now making Polygon Object
1578  std::pair<float,float> tmpvertex;
1579  //make Polygon Object as in mac/PolyOverlap.cc
1580  std::vector<std::pair<float,float> > vertices;
1581  for (unsigned int i=0; i<container_polygon.size(); i++){
1582  tmpvertex = std::make_pair( container_polygon.at(i)->w,
1583  container_polygon.at(i)->t );
1584  vertices.push_back( tmpvertex );
1585  }
1586  fParams.PolyObject = Polygon2D( vertices );
1587  }
1588 
1589  fTimeRecord_ProcName.push_back("FillPolygon");
1590  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1591  }
util::GeometryUtilities * fGSer
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 SelectPolygonHitList(const std::vector< util::PxHit > &hitlist, std::vector< const util::PxHit * > &hitlistlocal)
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 278 of file ClusterParamsAlg.cxx.

References lar::util::StatCollector< T, W >::add(), lar::util::StatCollector< T, W >::Average(), cluster::cluster_params::charge_wgt_x, cluster::cluster_params::charge_wgt_y, cluster::cluster_params::eigenvalue_principal, cluster::cluster_params::eigenvalue_secondary, fFinishedGetAverages, fHitVector, fParams, fTimeRecord_ProcName, fTimeRecord_ProcTime, cluster::cluster_params::mean_ADC, cluster::cluster_params::mean_charge, cluster::cluster_params::mean_x, cluster::cluster_params::mean_y, cluster::cluster_params::multi_hit_wires, cluster::cluster_params::N_Hits, cluster::cluster_params::N_Wires, lar::util::StatCollector< T, W >::RMS(), cluster::cluster_params::rms_ADC, cluster::cluster_params::rms_charge, lar::util::StatCollector< T, W >::Sum(), cluster::cluster_params::sum_ADC, and cluster::cluster_params::sum_charge.

Referenced by FillParams(), GetParams(), GetRoughAxis(), cluster::StandardClusterParamsAlg::Integral(), cluster::StandardClusterParamsAlg::IntegralStdDev(), cluster::StandardClusterParamsAlg::MultipleHitDensity(), MultipleHitDensity(), MultipleHitWires(), cluster::StandardClusterParamsAlg::NHits(), cluster::StandardClusterParamsAlg::SummedADC(), and cluster::StandardClusterParamsAlg::SummedADCStdDev().

278  {
279  if(!override) { //Override being set, we skip all this logic.
280  //OK, no override. Stop if we're already finshed.
281  if (fFinishedGetAverages) return;
282  }
283 
284  TStopwatch localWatch;
285  localWatch.Start();
286 
287  TPrincipal fPrincipal(2,"D");
288 
289  fParams.N_Hits = fHitVector.size();
290 
291  std::map<double, int> wireMap;
292 
293  lar::util::StatCollector<double> charge, sumADC;
294 
295  int uniquewires = 0;
296  int multi_hit_wires = 0;
297  for(auto& hit : fHitVector){
298  // std::cout << "This hit has charge " << hit.charge << "\n";
299 
300  double data[2];
301  data[0] = hit.w;
302  data[1] = hit.t;
303  fPrincipal.AddRow(data);
304  fParams.charge_wgt_x += hit.w * hit.charge;
305  fParams.charge_wgt_y += hit.t * hit.charge;
306  charge.add(hit.charge);
307  sumADC.add(hit.sumADC);
308 
309 
310  if (wireMap[hit.w] == 0) {
311  uniquewires ++;
312  }
313  if (wireMap[hit.w] == 1) {
314  multi_hit_wires ++;
315  }
316  wireMap[hit.w] ++;
317 
318 
319  }
320 
321  fParams.sum_charge = charge.Sum();
322  fParams.mean_charge = charge.Average();
323  fParams.rms_charge = charge.RMS();
324 
325  fParams.sum_ADC = sumADC.Sum();
326  fParams.mean_ADC = sumADC.Average();
327  fParams.rms_ADC = sumADC.RMS();
328 
329 
330  fParams.N_Wires = uniquewires;
331  fParams.multi_hit_wires = multi_hit_wires;
332 
333  if(fPrincipal.GetMeanValues()->GetNrows()<2) {
334  throw cluster::CRUException();
335  return;
336  }
337 
338  fParams.mean_x = (* fPrincipal.GetMeanValues())[0];
339  fParams.mean_y = (* fPrincipal.GetMeanValues())[1];
341 
342  if (fParams.sum_charge != 0.) {
345  }
346  else { // "SNAFU"; use the mean
349  }
350 
351  fPrincipal.MakePrincipals();
352 
353  fParams.eigenvalue_principal = (* fPrincipal.GetEigenValues() )[0];
354  fParams.eigenvalue_secondary = (* fPrincipal.GetEigenValues() )[1];
355 
356  fFinishedGetAverages = true;
357  // Report();
358 
359  fTimeRecord_ProcName.push_back("GetAverages");
360  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
361  }
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 ( bool  override_ = false)

Calculates the following variables: start_charge end_charge

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

Definition at line 1680 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::end_charge, EndCharge(), fFinishedGetEndCharges, fParams, fTimeRecord_ProcName, fTimeRecord_ProcTime, LinearIntegral(), cluster::cluster_params::start_charge, and StartCharge().

Referenced by FillParams(), and GetParams().

1680  {
1681  if(!override_) { //Override being set, we skip all this logic.
1682  //OK, no override. Stop if we're already finshed.
1683  if (fFinishedGetEndCharges) return;
1684  }
1685 
1686  TStopwatch localWatch;
1687  localWatch.Start();
1688 
1691 
1692  fFinishedGetEndCharges = true;
1693  // Report();
1694 
1695  fTimeRecord_ProcName.push_back("GetEndCharges");
1696  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1697  } // ClusterParamsAlg::GetEndCharges()
std::vector< std::string > fTimeRecord_ProcName
double start_charge
Charge at the start of the cluster.
Definition: ClusterParams.h:45
double EndCharge(float length=1., unsigned int nbins=10)
Returns the expected charge at the end of the cluster.
double StartCharge(float length=1., unsigned int nbins=10)
Returns the expected charge at the beginning of the cluster.
cluster::cluster_params fParams
std::vector< double > fTimeRecord_ProcTime
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 125 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::closing_angle, cluster::cluster_params::closing_angle_charge_wgt, cluster::cluster_params::eigenvalue_principal, cluster::cluster_params::eigenvalue_secondary, fParams, cluster::cluster_params::hit_density_1D, cluster::cluster_params::length, cluster::cluster_params::mean_charge, cluster::cluster_params::modified_hit_density, cluster::cluster_params::multi_hit_wires, cluster::cluster_params::N_Hits, cluster::cluster_params::N_Hits_HC, cluster::cluster_params::N_Wires, cluster::cluster_params::opening_angle, cluster::cluster_params::opening_angle_charge_wgt, PI, cluster::cluster_params::RMS_charge, cluster::cluster_params::sum_charge, and cluster::cluster_params::width.

Referenced by PrintFANNVector(), SetVerbose(), and TrackShowerSeparation().

125  {
126  unsigned int length = 13;
127  if (data.size() != length)
128  data.resize(length);
129  data[0] = fParams.opening_angle / PI;
130  data[1] = fParams.opening_angle_charge_wgt / PI;
131  data[2] = fParams.closing_angle / PI;
132  data[3] = fParams.closing_angle_charge_wgt / PI;
133  data[4] = fParams.eigenvalue_principal;
134  data[5] = fParams.eigenvalue_secondary;
135  data[6] = fParams.width / fParams.length;
138  data[9] = fParams.N_Hits_HC / fParams.N_Hits;
139  data[10] = fParams.modified_hit_density;
140  data[11] = fParams.RMS_charge / fParams.mean_charge;
141  data[12] = log(fParams.sum_charge / fParams.length);
142  return;
143  }
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
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
#define PI
double closing_angle
Width of angular distubtion wrt endpoint.
Definition: ClusterParams.h:43
void cluster::ClusterParamsAlg::GetFinalSlope ( 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 1130 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::angle_2d, cd(), Draw(), fChargeProfile, fChargeProfileNew, fCoarseChargeProfile, fCoarseNbins, fFinishedGetFinalSlope, fFinishedRefineStartPoints, fGSer, fHitVector, fParams, fProfileNbins, fProjectedLength, fRough2DIntercept, fRough2DSlope, fTimeRecord_ProcName, fTimeRecord_ProcTime, util::GeometryUtilities::Get2Dangle(), util::GeometryUtilities::Get2DDistance(), util::GeometryUtilities::GetPointOnLine(), cluster::cluster_params::hit_density_1D, cluster::cluster_params::length, cluster::cluster_params::mean_charge, cluster::cluster_params::modified_hit_density, cluster::cluster_params::modmeancharge, PI, RefineStartPoints(), SetLineColor(), cluster::cluster_params::start_point, Update(), and verbose.

Referenced by FillParams(), and GetParams().

1130  {
1131  if(!override) { //Override being set, we skip all this logic.
1132  //OK, no override. Stop if we're already finshed.
1133  if (fFinishedGetFinalSlope) return;
1134  //Try to run the previous function if not yet done.
1136  } else {
1137  //Try to run the previous function if not yet done.
1139  }
1140 
1141  TStopwatch localWatch;
1142  localWatch.Start();
1150  //double xangle=Get2DAngleForHit( wstn,tstn, hitlist);
1151  const int NBINS=720;
1152  std::vector< int > fh_omega_single(NBINS,0); //720,-180., 180.
1153 
1154  // fParams.start_point
1155 
1156  double current_maximum=0;
1157  double curr_max_bin=-1;
1158 
1159  for( auto hit : fHitVector){
1160 
1161  double omx=fGSer->Get2Dangle((util::PxPoint *)&hit,&fParams.start_point); // in rad and assuming cm/cm space.
1162  int nbin=(omx+TMath::Pi())*(NBINS-1)/(2*TMath::Pi());
1163  if(nbin >= NBINS) nbin=NBINS-1;
1164  if(nbin < 0) nbin=0;
1165  fh_omega_single[nbin]+=hit.charge;
1166 
1167  if( fh_omega_single[nbin]>current_maximum )
1168  {
1169  current_maximum= fh_omega_single[nbin];
1170  curr_max_bin=nbin;
1171  }
1172 
1173 
1174  }
1175 
1176 
1177  fParams.angle_2d=(curr_max_bin/720*(2*TMath::Pi()))-TMath::Pi();
1178  fParams.angle_2d*=180/PI;
1179  if(verbose) std::cout << " Final 2D angle: " << fParams.angle_2d << " degrees " << std::endl;
1180 
1181  double mod_angle=(fabs(fParams.angle_2d)<=90) ? fabs(fParams.angle_2d) : 180 - fabs(fParams.angle_2d); //want to transfer angle to radians and from 0 to 90.
1182 
1183  if(cos(mod_angle*PI/180.))
1184  { //values here based on fit of hit_density_1D vs. mod_angle defined as above (in ArgoNeuT).
1185  if(mod_angle<=75.)
1186  fParams.modified_hit_density=fParams.hit_density_1D/(2.45*cos(mod_angle*PI/180.)+0.2595);
1187  else
1188  fParams.modified_hit_density=fParams.hit_density_1D/(1.036*mod_angle*PI/180.-0.2561);
1189 
1190  //calculate modified mean_charge:
1191  fParams.modmeancharge = fParams.mean_charge/(1264 - 780*cos(mod_angle*PI/180.));
1192 
1193  }
1194  else
1196 
1198  // do not use for now.
1199  bool drawProfileHisto=false;
1200  if (drawProfileHisto)
1201  {
1202 
1204  double corr_factor=1;
1205  if(cos(mod_angle*PI/180.))
1206  { //values here based on fit of hit_density_1D vs. mod_angle defined as above (in ArgoNeuT).
1207 
1208  if(mod_angle<=75.)
1209  corr_factor*=(2.45*cos(mod_angle*PI/180.)+0.2595);
1210  else
1211  corr_factor*=(1.036*mod_angle*PI/180.-0.2561);
1212  }
1213 
1214  fProfileNbins=(int)(fProfileNbins/2*corr_factor + 0.5); // +0.5 to round to nearest sensible value
1215  //if(fProfileNbins<10) fProfileNbins=10;
1216 
1217  if(verbose) std::cout << " number of final profile bins " << fProfileNbins <<std::endl;
1218 
1219  fChargeProfile.clear();
1220  fCoarseChargeProfile.clear();
1221  fChargeProfile.resize(fProfileNbins,0);
1223 
1224 
1225  fChargeProfileNew.clear();
1226  // fDiffChargeProfile.clear();
1227  fChargeProfileNew.resize(fProfileNbins,0);
1228  // fDiffChargeProfile.resize(fProfileNbins,0);
1229 
1230 
1231  util::PxPoint BeginOnlinePoint;
1232 
1234 
1235  for( auto hit : fHitVector){
1236 
1237  util::PxPoint OnlinePoint;
1238  // get coordinates of point on axis.
1239  // std::cout << BeginOnlinePoint << std::endl;
1240  //std::cout << &OnlinePoint << std::endl;
1241  fGSer->GetPointOnLine(fRough2DSlope,&BeginOnlinePoint,&hit,OnlinePoint);
1242 
1243  double linedist=fGSer->Get2DDistance(&OnlinePoint,&BeginOnlinePoint);
1244  // double ortdist=fGSer->Get2DDistance(&OnlinePoint,&hit);
1245 
1246 
1247  int fine_bin =(int)(linedist/fProjectedLength*(fProfileNbins-1));
1248 
1249 
1250  if(fine_bin<fProfileNbins) //only fill if bin number is in range
1251  {
1252  fChargeProfileNew.at(fine_bin)+=hit.charge;
1253  }
1254  }
1255 
1256 
1257 
1258 
1259  TH1F * charge_histo =
1260  new TH1F("charge dist","charge dist",fProfileNbins,0,fProfileNbins);
1261 
1262  TH1F * charge_diff =
1263  new TH1F("charge diff","charge diff",fProfileNbins,0,fProfileNbins);
1264 
1265  for(int ix=0;ix<fProfileNbins-1;ix++)
1266  {
1267  charge_histo->SetBinContent(ix,fChargeProfileNew[ix]);
1268 
1269  if(ix>2 && ix < fProfileNbins-3)
1270  {
1271  double diff=( 1./12.*fChargeProfileNew[ix-2]
1272  - 2./3.*fChargeProfileNew [ix-1]
1273  + 2./3.*fChargeProfileNew [ix+1]
1274  - 1./12.*fChargeProfileNew[ix+2] )
1275  /(double)fProfileNbins;
1276  charge_diff->SetBinContent(ix,diff);
1277  }
1278  }
1279 
1280  TCanvas * chCanv = new TCanvas("chCanv","chCanv", 600, 600);
1281  chCanv -> cd();
1282  charge_histo -> SetLineColor(3);
1283  charge_histo -> Draw("");
1284  charge_diff -> SetLineColor(2);
1285  charge_diff -> Draw("same");
1286  chCanv -> Update();
1287 
1288 
1289 
1290 
1291  }
1292 
1294 
1295  fTimeRecord_ProcName.push_back("GetFinalSlope");
1296  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1297 
1298  fFinishedGetFinalSlope = true;
1299  return;
1300  }
util::GeometryUtilities * fGSer
std::vector< double > fChargeProfileNew
std::vector< std::string > fTimeRecord_ProcName
Double_t Get2DDistance(Double_t wire1, Double_t time1, Double_t wire2, Double_t time2) const
void RefineStartPoints(bool override=false)
Double_t Get2Dangle(Double_t deltawire, Double_t deltatime) const
std::vector< util::PxHit > fHitVector
gr2 SetLineColor(2)
c1 Update()
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.
Int_t GetPointOnLine(Double_t slope, Double_t intercept, Double_t wire1, Double_t time1, Double_t &wireout, Double_t &timeout) const
std::vector< double > fCoarseChargeProfile
double mean_charge
Mean (average) charge of hits in ADC.
Definition: ClusterParams.h:28
hist1 Draw("HIST")
c1 cd(1)
#define PI
std::vector< double > fChargeProfile
const std::vector<util::PxHit>& cluster::ClusterParamsAlg::GetHitVector ( ) const
inline

Definition at line 271 of file ClusterParamsAlg.h.

References fHitVector.

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

Definition at line 270 of file ClusterParamsAlg.h.

References fHitVector.

Referenced by cluster::StandardClusterParamsAlg::NInputHits(), and cmtool::CPriorityAlgoBase::Priority().

270 {return fHitVector.size();}
std::vector< util::PxHit > fHitVector
void cluster::ClusterParamsAlg::GetOpeningAngle ( )

Referenced by setNeuralNetPath().

void cluster::ClusterParamsAlg::GetProfileInfo ( 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 437 of file ClusterParamsAlg.cxx.

References AddEntry(), cd(), Draw(), cluster::cluster_params::end_point, fBeginIntercept, fChargeCutoffThreshold, fChargeProfile, fCoarseChargeProfile, fCoarseNbins, fEndIntercept, fFinishedGetProfileInfo, fFinishedGetRoughAxis, fGSer, fHitVector, Fill(), fInterHigh_side, fInterLow_side, fParams, fPlane, fProfileIntegralBackward, fProfileIntegralForward, fProfileMaximumBin, fProfileNbins, fProjectedLength, fRough2DIntercept, fRough2DSlope, fRoughBeginPoint, fRoughEndPoint, fTimeRecord_ProcName, fTimeRecord_ProcTime, util::GeometryUtilities::Get2DDistance(), util::GeometryUtilities::GetPointOnLine(), util::GeometryUtilities::GetPointOnLineWSlopes(), GetRoughAxis(), leg, cluster::cluster_params::N_Hits, util::PxPoint::plane, Scale(), SetLineColor(), cluster::cluster_params::start_point, cluster::cluster_params::sum_charge, util::PxPoint::t, Update(), verbose, util::PxPoint::w, weight, and cluster::cluster_params::width.

Referenced by EndCharge(), FillParams(), GetParams(), IntegrateFitCharge(), RefineStartPointAndDirection(), and cluster::StandardClusterParamsAlg::Width().

437  {
438  if(!override) { //Override being set, we skip all this logic.
439  //OK, no override. Stop if we're already finshed.
440  if (fFinishedGetProfileInfo) return;
441  //Try to run the previous function if not yet done.
443  } else {
445  }
446 
447  TStopwatch localWatch;
448  localWatch.Start();
449 
450  bool drawOrtHistos = false;
451 
452  //these variables need to be initialized to other values?
453  if(fRough2DSlope==-999.999 || fRough2DIntercept==-999.999 )
454  GetRoughAxis(true);
455 
456  //get slope of lines orthogonal to those found crossing the shower.
457  double inv_2d_slope=0;
458  if(fabs(fRough2DSlope)>0.001)
459  inv_2d_slope=-1./fRough2DSlope;
460  else
461  inv_2d_slope=-999999.;
462  // std::cout << " N_Hits is " << fParams.N_Hits << std::endl;
463  // std::cout << " inv_2d_slope is " << inv_2d_slope << std::endl;
464  // std::cout << " fRough2DSlope is " << fRough2DSlope << std::endl;
465  // std::cout << " conversions " << fWire2Cm.at(fPlane)
466  // << " " << fTime2Cm << std::endl;
467 
468  double InterHigh=-999999;
469  double InterLow=999999;
470  double WireHigh=-999999;
471  double WireLow=999999;
472  fInterHigh_side=-999999;
473  fInterLow_side=999999;
474  double min_ortdist(999), max_ortdist(-999); // needed to calculate width
475  //loop over all hits. Create coarse and fine profiles
476  // of the charge weight to help refine the start/end
477  // points and have a first guess of direction
478 
479  //std::cout << " inv_2d_slope " << inv_2d_slope << std::endl;
480 
481  for(auto const &hit : fHitVector)
482  {
483 
484  //calculate intercepts along
485  double intercept=hit.t - inv_2d_slope * (double)(hit.w);
486  double side_intercept=hit.t - fRough2DSlope * (double)(hit.w);
487 
488  util::PxPoint OnlinePoint;
489  // get coordinates of point on axis.
491 
492  // std::cout << "normal point " << hit.w << ","<<hit.t << " online point " << OnlinePoint.w << "," << OnlinePoint.t << std::endl;
493 
494  double ortdist=fGSer->Get2DDistance(&OnlinePoint,&hit);
495  if (ortdist < min_ortdist) min_ortdist = ortdist;
496  if (ortdist > max_ortdist) max_ortdist = ortdist;
497 
498  if(inv_2d_slope!=-999999) //almost all cases
499  {
500  if(intercept > InterHigh ){
501  InterHigh=intercept;
502  }
503 
504  if(intercept < InterLow ){
505  InterLow=intercept;
506  }
507  }
508  else //slope is practically horizontal. Care only about wires.
509  {
510  if(hit.w > WireHigh ){
511  WireHigh=hit.w;
512  }
513 
514  if(hit.w < WireLow ){
515  WireLow=hit.w;
516  }
517  }
518 
519  if(side_intercept > fInterHigh_side ){
520  fInterHigh_side=side_intercept;
521  }
522 
523  if(side_intercept < fInterLow_side ){
524  fInterLow_side=side_intercept;
525  }
526 
527 
528  }
529  // end of first HitIter loop, at this point we should
530  // have the extreme intercepts
531 
533  // Second loop. Fill profiles.
535 
536  // get projected points at the beginning and end of the axis.
537 
538  util::PxPoint HighOnlinePoint, LowOnlinePoint,BeginOnlinePoint, EndOnlinePoint;
539 
540  if(inv_2d_slope!=-999999) //almost all cases
541  {
543  InterHigh,HighOnlinePoint);
545  InterLow,LowOnlinePoint);
546  }
547  else //need better treatment of horizontal events.
548  {
549  util::PxPoint ptemphigh(fPlane,WireHigh,(fInterHigh_side+fInterLow_side)/2);
550  util::PxPoint ptemplow(fPlane,WireLow,(fInterHigh_side+fInterLow_side)/2);
551  fGSer->GetPointOnLine(fRough2DSlope,fRough2DIntercept,&ptemphigh,HighOnlinePoint);
552  fGSer->GetPointOnLine(fRough2DSlope,fRough2DIntercept,&ptemplow,LowOnlinePoint);
553  }
554  if(verbose){
555 
556  // std::cout << " extreme intercepts: low: " << InterLow
557  // << " " << InterHigh << std::endl;
558  // std::cout << " extreme intercepts: side: " << fInterLow_side
559  // << " " << fInterHigh_side << std::endl;
560  // std::cout << "axis + intercept " << fRough2DSlope << " "
561  // << fRough2DIntercept << std::endl;
562  //
563  std::cout << " Low online point: " << LowOnlinePoint.w << ", " << LowOnlinePoint.t
564  << " High: " << HighOnlinePoint.w << ", " << HighOnlinePoint.t << std::endl;
565 
566  //std::cout << " max_min ortdist" << max_ortdist << " " << min_ortdist << std::endl;
567  //std::cout << " hit list size " << fHitVector.size() << std::endl;
568  //define BeginOnlinePoint as the one with lower wire number (for now), adjust intercepts accordingly
569  }
570  if(HighOnlinePoint.w >= LowOnlinePoint.w)
571  {
572  BeginOnlinePoint=LowOnlinePoint;
573  fBeginIntercept=InterLow;
574  EndOnlinePoint=HighOnlinePoint;
575  fEndIntercept=InterHigh;
576  }
577  else
578  {
579  BeginOnlinePoint=HighOnlinePoint;
580  fBeginIntercept=InterHigh;
581  EndOnlinePoint=LowOnlinePoint;
582  fEndIntercept=InterLow;
583  }
584 
585  fProjectedLength=fGSer->Get2DDistance(&BeginOnlinePoint,&EndOnlinePoint);
586 
587  // std::cout << " projected length " << fProjectedLength
588  // << " Begin Point " << BeginOnlinePoint.w << " "
589  // << BeginOnlinePoint.t << " End Point " << EndOnlinePoint.w << ","<< EndOnlinePoint.t << std::endl;
590 
592  fCoarseNbins=2;
593 
595  if(fProfileNbins<10) fProfileNbins=10;
596  //std::cout << " number of profile bins " << fProfileNbins <<std::endl;
597 
598  fChargeProfile.clear();
599  fCoarseChargeProfile.clear();
600  fChargeProfile.resize(fProfileNbins,0);
602 
603 
605  // Some fitting variables to make a histogram:
606 
607  // TODO this is nonsense for small clusters
608  std::vector<double> ort_profile;
609  const int NBINS=100;
610  ort_profile.resize(NBINS);
611 
612  std::vector<double> ort_dist_vect;
613  ort_dist_vect.reserve(fHitVector.size());
614 
615  double current_maximum=0;
616  for(auto& hit : fHitVector)
617  {
618 
619  util::PxPoint OnlinePoint;
620  // get coordinates of point on axis.
621  // std::cout << BeginOnlinePoint << std::endl;
622  //std::cout << &OnlinePoint << std::endl;
623  fGSer->GetPointOnLine(fRough2DSlope,&BeginOnlinePoint,&hit,OnlinePoint);
624  double ortdist=fGSer->Get2DDistance(&OnlinePoint,&hit);
625 
626  double linedist=fGSer->Get2DDistance(&OnlinePoint,&BeginOnlinePoint);
627  // double ortdist=fGSer->Get2DDistance(&OnlinePoint,&hit);
628  ort_dist_vect.push_back(ortdist);
629  int ortbin;
630  if (ortdist == 0)
631  ortbin = 0;
632  else if (max_ortdist == min_ortdist)
633  ortbin = 0;
634  else
635  ortbin =(int)(ortdist-min_ortdist)/(max_ortdist-min_ortdist)*(NBINS-1);
636 
637  ort_profile.at(ortbin)+=hit.charge;
638  //if (ortdist < min_ortdist) min_ortdist = ortdist;
639  //if (ortdist > max_ortdist) max_ortdist = ortdist;
640 
642  //calculate the weight along the axis, this formula is based on rough guessology.
643  // there is no physics motivation behind the particular numbers, A.S.
644  // \todo: switch to something that is motivated and easier to
645  // spell than guessology. C.A.
647  double weight= (ortdist<1.) ? 10 * (hit.charge) : 5 * (hit.charge) / ortdist;
648 
649  int fine_bin = fProjectedLength?
650  (int)(linedist/fProjectedLength*(fProfileNbins-1)): 0;
651  int coarse_bin= fProjectedLength?
652  (int)(linedist/fProjectedLength*(fCoarseNbins-1)): 0;
653  /*
654  std::cout << "linedist: " << linedist << std::endl;
655  std::cout << "fProjectedLength: " << fProjectedLength << std::endl;
656  std::cout << "fProfileNbins: " << fProfileNbins << std::endl;
657  std::cout << "fine_bin: " << fine_bin << std::endl;
658  std::cout << "coarse_bin: " << coarse_bin << std::endl;
659  */
660 
661  // std::cout << "length" << linedist << " fine_bin, coarse " << fine_bin << " " << coarse_bin << std::endl;
662 
663 
664 
665  if(fine_bin<fProfileNbins) //only fill if bin number is in range
666  {
667  fChargeProfile.at(fine_bin)+=weight;
668 
669  //find maximum bin on the fly:
670  if(fChargeProfile.at(fine_bin)>current_maximum
671  && fine_bin!=0 && fine_bin!=fProfileNbins-1)
672  {
673  current_maximum=fChargeProfile.at(fine_bin);
674  fProfileMaximumBin=fine_bin;
675  }
676  }
677 
678  if(coarse_bin<fCoarseNbins) //only fill if bin number is in range
679  fCoarseChargeProfile.at(coarse_bin)+=weight;
680 
681  } // end second loop on hits. Now should have filled profile vectors.
682 
683  if(verbose) std::cout << "end second loop " << std::endl;
684 
685  double integral=0;
686  int ix=0;
687  // int fbin=0;
688  for(ix=0;ix<NBINS;ix++)
689  {
690  integral+=ort_profile.at(ix);
691  if(integral>=0.95*fParams.sum_charge)
692  {
693  break;
694  }
695  }
696 
697  fParams.width=2*(double)ix/(double)(NBINS-1)*(double)(max_ortdist-min_ortdist); // multiply by two because ortdist is folding in both sides.
698 
699  if(verbose) std::cout << " after width " << std::endl;
700 
701 
702  if (drawOrtHistos){
703  TH1F * ortDistHist =
704  new TH1F("ortDistHist",
705  "Orthogonal Distance to axis;Distance (cm)",
706  100, min_ortdist, max_ortdist);
707  TH1F * ortDistHistCharge =
708  new TH1F("ortDistHistCharge",
709  "Orthogonal Distance to axis (Charge Weighted);Distance (cm)",
710  100, min_ortdist, max_ortdist);
711  TH1F * ortDistHistAndrzej=
712  new TH1F("ortDistHistAndrzej",
713  "Orthogonal Distance Andrzej weighting",
714  100, min_ortdist, max_ortdist);
715 
716  for (int index = 0; index < fParams.N_Hits; index++){
717  ortDistHist -> Fill(ort_dist_vect.at(index));
718  ortDistHistCharge -> Fill(ort_dist_vect.at(index), fHitVector.at(index).charge);
719  double weight= (ort_dist_vect.at(index)<1.) ?
720  10 * (fHitVector.at(index).charge) :
721  (5 * (fHitVector.at(index).charge)/ort_dist_vect.at(index));
722  ortDistHistAndrzej -> Fill(ort_dist_vect.at(index), weight);
723  }
724  ortDistHist -> Scale(1.0/ortDistHist->Integral());
725  ortDistHistCharge -> Scale(1.0/ortDistHistCharge->Integral());
726  ortDistHistAndrzej -> Scale(1.0/ortDistHistAndrzej->Integral());
727 
728  TCanvas * ortCanv = new TCanvas("ortCanv","ortCanv", 600, 600);
729  ortCanv -> cd();
730  ortDistHistAndrzej -> SetLineColor(3);
731  ortDistHistAndrzej -> Draw("");
732  ortDistHistCharge -> Draw("same");
733  ortDistHist -> SetLineColor(2);
734  ortDistHist -> Draw("same");
735 
736  TLegend * leg = new TLegend(.4,.6,.8,.9);
737  leg -> SetHeader("Charge distribution");
738  leg -> AddEntry(ortDistHist, "Unweighted Hits");
739  leg -> AddEntry(ortDistHistCharge, "Charge Weighted Hits");
740  leg -> AddEntry(ortDistHistAndrzej, "Charge Weighted by Guessology");
741  leg -> Draw();
742 
743  ortCanv -> Update();
744  }
745 
748 
749  //calculate the forward and backward integrals counting int the maximum bin.
750 
751 
752  for(int ibin=0;ibin<fProfileNbins;ibin++)
753  {
754  if(ibin<=fProfileMaximumBin)
756 
757  if(ibin>=fProfileMaximumBin)
759  }
760 
761  // now, we have the forward and backward integral so start
762  // stepping forward and backward to find the trunk of the
763  // shower. This is done making sure that the running
764  // integral until less than 1% is left and the bin is above
765  // a set theshold many empty bins.
766 
767 
768 
769  //forward loop
770  double running_integral=fProfileIntegralForward;
771  int startbin,endbin;
772  for(startbin=fProfileMaximumBin; startbin>1 && startbin < (int)(fChargeProfile.size());startbin--)
773  {
774  running_integral-=fChargeProfile.at(startbin);
775  if( fChargeProfile.at(startbin)<fChargeCutoffThreshold.at(fPlane) && running_integral/fProfileIntegralForward<0.01 )
776  break;
777  else if(fChargeProfile.at(startbin)<fChargeCutoffThreshold.at(fPlane)
778  && startbin-1>0 && fChargeProfile.at(startbin-1)<fChargeCutoffThreshold.at(fPlane)
779  && startbin-2>0 && fChargeProfile.at(startbin-2)<fChargeCutoffThreshold.at(fPlane))
780  break;
781  }
782 
783  //backward loop
784  running_integral=fProfileIntegralBackward;
785  for(endbin=fProfileMaximumBin; endbin>0 && endbin<fProfileNbins-1; endbin++)
786  {
787  running_integral-=fChargeProfile.at(endbin);
788  if( fChargeProfile.at(endbin)<fChargeCutoffThreshold.at(fPlane) && running_integral/fProfileIntegralBackward<0.01 )
789  break;
790  else if(fChargeProfile.at(endbin)<fChargeCutoffThreshold.at(fPlane)
791  && endbin+1<fProfileNbins-1 && endbin+2<fProfileNbins-1
792  && fChargeProfile.at(endbin+1)<fChargeCutoffThreshold.at(fPlane)
793  && fChargeProfile.at(endbin+2)<fChargeCutoffThreshold.at(fPlane))
794  break;
795  }
796 
797  // now have profile start and endpoints. Now translate to wire/time.
798  // Will use wire/time that are on the rough axis.
799  // fProjectedLength is the length from fInterLow to interhigh a
800  // long the rough_2d_axis on bin distance is:
801  // util::PxPoint OnlinePoint;
802 
803  if(inv_2d_slope!=-999999) //almost all cases
804  {
805  double ort_intercept_begin=fBeginIntercept+(fEndIntercept-fBeginIntercept)/fProfileNbins*startbin;
806  //std::cout << " ort_intercept_begin: " << ort_intercept_begin << std::endl;
809  ort_intercept_begin,
811 
812  double ort_intercept_end=fBeginIntercept+(fEndIntercept-fBeginIntercept)/fProfileNbins*endbin;
814 
815  //std::cout << " ort_intercept_end: " << ort_intercept_end << std::endl;
818  ort_intercept_end,
821  }
822  else
823  {
824  double wire_begin=WireLow+(WireHigh-WireLow)/fProfileNbins*startbin;
825  //std::cout << " wire_begin: " << wire_begin << std::endl;
826 
827  util::PxPoint ptemphigh(fPlane,wire_begin,(fInterHigh_side+fInterLow_side)/2);
830 
831  double wire_end=WireLow+(WireHigh-WireLow)/fProfileNbins*endbin;
832  //std::cout << " wire_end: " << wire_end << std::endl;
833 
834  util::PxPoint ptemplow(fPlane,wire_end,(fInterHigh_side+fInterLow_side)/2);
837 
838  }
839 
840  if(verbose) std::cout << " rough start points " << fRoughBeginPoint.w << ", " << fRoughBeginPoint.t << " end: " << fRoughEndPoint.w << " " << fRoughEndPoint.t << std::endl;
841 
842  // ok. now have fRoughBeginPoint and fRoughEndPoint. No decision about direction has been made yet.
845  // fRoughEndPoint
846  // fRoughEndPoint
847  // and use them to get the axis
848 
849 
851  // Report();
852 
853  fTimeRecord_ProcName.push_back("GetProfileInfo");
854  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
855 
856  return;
857  }
util::GeometryUtilities * fGSer
void GetRoughAxis(bool override=false)
std::vector< std::string > fTimeRecord_ProcName
Double_t Get2DDistance(Double_t wire1, Double_t time1, Double_t wire2, Double_t time2) const
TNtupleSim Fill(f1, f2, f3, f4)
TLegend * leg
Definition: compare.C:67
std::vector< util::PxHit > fHitVector
gr2 SetLineColor(2)
c1 Update()
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
Int_t GetPointOnLine(Double_t slope, Double_t intercept, Double_t wire1, Double_t time1, Double_t &wireout, Double_t &timeout) const
leg AddEntry(gr1,"Reference data","lp")
double weight
Definition: plottest35.C:25
std::vector< double > fCoarseChargeProfile
util::PxPoint end_point
end point
Definition: ClusterParams.h:25
hist1 Draw("HIST")
Int_t GetPointOnLineWSlopes(Double_t slope, Double_t intercept, Double_t ort_intercept, Double_t &wireout, Double_t &timeout) const
c1 cd(1)
simulatedPeak Scale(1/simulationNormalisationFactor)
std::vector< double > fChargeProfile
unsigned int plane
Definition: PxUtils.h:12
std::vector< double > fChargeCutoffThreshold
void cluster::ClusterParamsAlg::GetRoughAxis ( bool  override = false)

Calculates the following variables: verticalness fRough2DSlope fRough2DIntercept

Parameters
override[description]

Definition at line 365 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::charge_wgt_x, cluster::cluster_params::charge_wgt_y, cluster::cluster_params::cluster_angle_2d, fFinishedGetAverages, fFinishedGetRoughAxis, fHitVector, fParams, fRough2DIntercept, fRough2DSlope, fTimeRecord_ProcName, fTimeRecord_ProcTime, GetAverages(), geo::vect::isfinite(), cluster::cluster_params::mean_charge, cluster::cluster_params::mean_x, cluster::cluster_params::mean_y, cluster::cluster_params::N_Hits, cluster::cluster_params::N_Hits_HC, PI, cluster::cluster_params::RMS_charge, cluster::cluster_params::rms_x, and cluster::cluster_params::rms_y.

Referenced by FillParams(), GetParams(), GetProfileInfo(), and cluster::StandardClusterParamsAlg::StartAngle().

365  {
366  if(!override) { //Override being set, we skip all this logic.
367  //OK, no override. Stop if we're already finshed.
368  if (fFinishedGetRoughAxis) return;
369  //Try to run the previous function if not yet done.
370  if (!fFinishedGetAverages) GetAverages(true);
371  } else {
372  //Try to run the previous function if not yet done.
373  if (!fFinishedGetAverages) GetAverages(true);
374  }
375 
376  TStopwatch localWatch;
377  localWatch.Start();
378 
379  double rmsx = 0.0;
380  double rmsy = 0.0;
381  double rmsq = 0.0;
382  //using the charge weighted coordinates find the axis from slope
383  double ncw=0;
384  double sumtime=0;//from sum averages
385  double sumwire=0;//from sum averages
386  double sumwiretime=0;//sum over (wire*time)
387  double sumwirewire=0;//sum over (wire*wire)
388  //next loop over all hits again
389 
390  for (auto& hit : fHitVector){
391  // First, abuse this loop to calculate rms in x and y
392  rmsx += pow(fParams.mean_x - hit.w, 2)/fParams.N_Hits;
393  rmsy += pow(fParams.mean_y - hit.t, 2)/fParams.N_Hits;
394  rmsq += pow(fParams.mean_charge - hit.charge, 2)/fParams.N_Hits;
395  //if charge is above avg_charge
396  // std::cout << "This hit has charge " << hit . charge << "\n";
397 
398  if(hit.charge > fParams.mean_charge){
399  ncw+=1;
400  sumwire+=hit.w;
401  sumtime+=hit.t;
402  sumwiretime+=hit.w * hit.t;
403  sumwirewire+=pow(hit.w,2);
404  }//for high charge
405  }//For hh loop
406 
407  fParams.rms_x = sqrt(rmsx);
408  fParams.rms_y = sqrt(rmsy);
409  fParams.RMS_charge = sqrt(rmsq);
410 
411  fParams.N_Hits_HC = ncw;
412  //Looking for the slope and intercept of the line above avg_charge hits
413 
414  // std::cout << " ncw, and parentheses " << ncw << " " << (ncw*sumwiretime- sumwire*sumtime) << " " <<(ncw*sumwirewire-sumwire*sumwire) << std::endl;
415  if((ncw*sumwirewire-sumwire*sumwire) > 0.00001)
416  fRough2DSlope= (ncw*sumwiretime- sumwire*sumtime)/(ncw*sumwirewire-sumwire*sumwire);//slope for cw
417  else
418  fRough2DSlope=999;
421  fParams.charge_wgt_y -fRough2DSlope*(fParams.charge_wgt_x): //intercept for cw
422  0.;
423 
424  //Getthe 2D_angle
426 
427 
428  fFinishedGetRoughAxis = true;
429  // Report();
430 
431  fTimeRecord_ProcName.push_back("GetRoughAxis");
432  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
433  return;
434  }
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
std::vector< util::PxHit > fHitVector
double mean_y
Mean of hits along y, peaks only.
Definition: ClusterParams.h:34
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
#define PI
void cluster::ClusterParamsAlg::Initialize ( )

Definition at line 181 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::Clear(), fChargeCutoffThreshold, fFinishedGetAverages, fFinishedGetEndCharges, fFinishedGetFinalSlope, fFinishedGetProfileInfo, fFinishedGetRoughAxis, fFinishedRefineDirection, fFinishedRefineStartPointAndDirection, fFinishedRefineStartPoints, fFinishedTrackShowerSep, fGSer, fHitVector, fParams, fProfileIntegralBackward, fProfileIntegralForward, fProfileMaximumBin, fRough2DIntercept, fRough2DSlope, fRoughBeginPoint, fRoughEndPoint, fTimeRecord_ProcName, fTimeRecord_ProcTime, util::GeometryUtilities::GetME(), util::PxPoint::t, and util::PxPoint::w.

Referenced by cluster::StandardClusterParamsAlg::Clear(), ClusterParamsAlg(), SetHits(), and ~ClusterParamsAlg().

182  {
183 
184  fTimeRecord_ProcName.clear();
185  fTimeRecord_ProcTime.clear();
186 
187  TStopwatch localWatch;
188  localWatch.Start();
189 
190  // Set pointer attributes
192 
193  //--- Initilize attributes values ---//
194  fFinishedGetAverages = false;
195  fFinishedGetRoughAxis = false;
196  fFinishedGetProfileInfo = false;
198  fFinishedRefineDirection = false;
199  fFinishedGetFinalSlope = false;
201  fFinishedTrackShowerSep = false;
202  fFinishedGetEndCharges = false;
203 
204  fRough2DSlope=-999.999; // slope
205  fRough2DIntercept=-999.999; // slope
206 
207  fRoughBeginPoint.w=-999.999;
208  fRoughEndPoint.w=-999.999;
209 
210  fRoughBeginPoint.t=-999.999;
211  fRoughEndPoint.t=-999.999;
212 
213  fProfileIntegralForward=999.999;
214  fProfileIntegralBackward=999.999;
215  fProfileMaximumBin=-999;
216 
217  fChargeCutoffThreshold.clear();
218  fChargeCutoffThreshold.reserve(3);
219  //fChargeCutoffThreshold.resize(3,0);
220  //fChargeCutoffThreshold.at(0)=200;
221  //fChargeCutoffThreshold.at(1)=400;
222  fChargeCutoffThreshold.push_back(500);
223  fChargeCutoffThreshold.push_back(500);
224  fChargeCutoffThreshold.push_back(1000);
225 
226  fHitVector.clear();
227 
228  fParams.Clear();
229 
230  // Initialize the neural network:
231  // enableFANN = false;
232  fTimeRecord_ProcName.push_back("Initialize");
233  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
234  }
util::GeometryUtilities * fGSer
static const GeometryUtilities * GetME()
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 ( 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 1709 of file ClusterParamsAlg.cxx.

References lar::util::details::SimplePolyFitterDataBase< T, D >::add(), fChargeProfile, lar::util::details::SimplePolyFitterBase< T, 1U >::FitParameters(), fProjectedLength, GetProfileInfo(), LinearIntegral(), min, and StartCharge().

Referenced by EndCharge(), and StartCharge().

1712  {
1713  // first compute the information on the charge profile
1714  GetProfileInfo();
1715 
1716  // first things first
1717  if (fit_first_bin > fit_end_bin) std::swap(fit_first_bin, fit_end_bin);
1718 
1719  // how many bins can we use?
1720  const unsigned int nbins = std::min
1721  (fit_end_bin - fit_first_bin, (unsigned int) fChargeProfile.size());
1722  if (nbins == 0) return 0;
1723 
1724  // move the specified range within reason
1725  if (fit_end_bin > fChargeProfile.size()) {
1726  fit_end_bin = fChargeProfile.size();
1727  fit_first_bin = fit_end_bin - nbins;
1728  }
1729 
1730  // if we have to use only one bin, so be it
1731  if (nbins < 2) return fChargeProfile[fit_first_bin];
1732 
1733  // fit the charge profile vs. bin number;
1734  // we assume that the first bin (#0) starts from the very beginning of the
1735  // cluster, that is at axis coordinate 0.
1737  for (unsigned int iBin = fit_first_bin; iBin < fit_end_bin; ++iBin) {
1738  // should we use a Poisson error instead of no error?
1739  fitter.add((double) iBin, fChargeProfile[iBin]);
1740  } // for
1741 
1742  // get the linear fit parameters; [0] intercept [1] slope
1744  try {
1745  fit = fitter.FitParameters();
1746  }
1747  catch (std::range_error) { // oops...
1748  // this is actually unexpected, since the only reason for the fit to fail
1749  // would be a singular fit matrix, that should be pretty much impossible
1750  // given that the bin coordinates are well behaved and there are at least
1751  // two of them
1752  std::cerr << "IntegrateFitCharge(): linear fit failed!" << std::endl;
1753  return 0.;
1754  }
1755 
1756  // now determine the bins corresponding to the length to integrate;
1757  // note that length can be pathologic (0, negative...); not our problem!
1758  const double length_to_bin
1759  = (double(fChargeProfile.size() - 1) / fProjectedLength);
1760  const double from_bin = from_length * length_to_bin,
1761  to_bin = to_length * length_to_bin;
1762 
1763  // we want the integral between from_bin and to_bin now:
1764  return LinearIntegral(fit[1] /* m */, fit[0] /* q */, from_bin, to_bin);
1765  } // ClusterParamsAlg::IntegrateFitCharge()
void GetProfileInfo(bool override=false)
Performs a linear regression of data.
Definition: SimpleFits.h:849
std::array< Data_t, NParams > FitParameters_t
type of set of fit parameters
Definition: SimpleFits.h:866
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
virtual FitParameters_t FitParameters() const override
Computes and returns all the parameters of the fit result.
Int_t min
Definition: plot.C:26
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].
std::vector< double > fChargeProfile
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 1702 of file ClusterParamsAlg.cxx.

References geo::sqr().

Referenced by GetEndCharges(), and IntegrateFitCharge().

1703  {
1704  return m / 2. * (sqr(x2) - sqr(x1)) + q * (x2 - x1);
1705  } // ClusterParamsAlg::LinearIntegral()
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqr(T v)
Float_t x2[n_points_geant4]
Definition: compare.C:26
size_t cluster::ClusterParamsAlg::MinNHits ( ) const
inline

Definition at line 63 of file ClusterParamsAlg.h.

References fMinNHits, and SetHits().

63 { return fMinNHits; }
size_t fMinNHits
Cut value for # hits: below this value clusters are not evaluated.
float cluster::ClusterParamsAlg::MultipleHitDensity ( )

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 1827 of file ClusterParamsAlg.cxx.

References fHitVector, fParams, GetAverages(), cluster::cluster_params::length, cluster::cluster_params::multi_hit_wires, and RefineStartPoints().

Referenced by RoughIntercept().

1827  {
1828  if (fHitVector.size() < 2) return 0.0F;
1829 
1830  // compute all the averages
1831  GetAverages();
1832  RefineStartPoints(); // fParams.length
1833 
1834  // return the relevant information
1835  return std::isnormal(fParams.length)?
1837  } // ClusterParamsAlg::MultipleHitDensity()
void RefineStartPoints(bool override=false)
std::vector< util::PxHit > fHitVector
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 1814 of file ClusterParamsAlg.cxx.

References fHitVector, fParams, GetAverages(), cluster::cluster_params::multi_hit_wires, and cluster::cluster_params::N_Wires.

Referenced by RoughIntercept().

1814  {
1815  if (fHitVector.size() < 2) return 0.0F;
1816 
1817  // compute all the averages
1818  GetAverages();
1819 
1820  // return the relevant information
1821  return std::isnormal(fParams.N_Wires)?
1823  } // 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 272 of file ClusterParamsAlg.h.

References fPlane, and SetPlane().

void cluster::ClusterParamsAlg::PrintFANNVector ( )

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

Returns
[description]

Definition at line 151 of file ClusterParamsAlg.cxx.

References GetFANNVector(), and verbose.

Referenced by SetVerbose().

151  {
152  std::vector<float> data;
153  GetFANNVector(data);
154  if(verbose){
155  std::cout << "Printing FANN input vector:\n"
156  << " Opening Angle (normalized) ... : " << data[0] << "\n"
157  << " Opening Angle charge weight .. : " << data[1] << "\n"
158  << " Closing Angle (normalized) ... : " << data[2] << "\n"
159  << " Closing Angle charge weight .. : " << data[3] << "\n"
160  << " Principal Eigenvalue ......... : " << data[4] << "\n"
161  << " Secondary Eigenvalue ......... : " << data[5] << "\n"
162  << " Width / Length ............... : " << data[6] << "\n"
163  << " Hit Density / M.H.D. ......... : " << data[7] << "\n"
164  << " Percent MultiHit Wires ....... : " << data[8] << "\n"
165  << " Percent High Charge Hits ..... : " << data[9] << "\n"
166  << " Modified Hit Density ......... : " << data[10] << "\n"
167  << " Charge RMS / Mean Charge ...... : " << data[11] << "\n"
168  << " log(Sum Charge / Length) ...... : " << data[12] << "\n";
169  }
170  return;
171  }
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 1308 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::closing_angle, cluster::cluster_params::closing_angle_charge_wgt, cluster::cluster_params::end_point, fCoarseChargeProfile, fFinishedRefineDirection, fFinishedRefineStartPoints, fHitVector, fParams, fQMinRefDir, fRoughBeginPoint, fRoughEndPoint, fTimeRecord_ProcName, fTimeRecord_ProcTime, cluster::cluster_params::mean_charge, cluster::cluster_params::N_Hits, cluster::cluster_params::N_Hits_HC, cluster::cluster_params::opening_angle, cluster::cluster_params::opening_angle_charge_wgt, PI, cluster::cluster_params::start_point, cluster::cluster_params::sum_charge, util::PxPoint::t, verbose, and util::PxPoint::w.

Referenced by cluster::StandardClusterParamsAlg::EndOpeningAngle(), GetParams(), RefineStartPointAndDirection(), and cluster::StandardClusterParamsAlg::StartOpeningAngle().

1308  {
1309  //
1310  // We don't use "override"? Should we remove? 05/01/14
1311  //
1312  if(!override) override = true;
1313 
1314  TStopwatch localWatch;
1315  localWatch.Start();
1316 
1317  // if(!override) { //Override being set, we skip all this logic.
1318  // //OK, no override. Stop if we're already finshed.
1319  // if (fFinishedRefineDirection) return;
1320  // //Try to run the previous function if not yet done.
1321  // if (!fFinishedGetProfileInfo) GetProfileInfo(true);
1322  // } else {
1323  // //Try to run the previous function if not yet done.
1324  // if (!fFinishedGetProfileInfo) GetProfileInfo(true);
1325  // }
1326 
1327  // double wire_2_cm = fGSer->WireToCm();
1328  // double time_2_cm = fGSer->TimeToCm();
1329 
1330  // Removing the conversion since these points are set above in cm-cm space
1331  // -- Corey
1332 
1333  util::PxPoint this_startPoint, this_endPoint;
1334 
1336  this_startPoint = fParams.start_point;
1337  this_endPoint = fParams.end_point;
1338  }
1339  else{
1340  this_startPoint = fRoughBeginPoint;
1341  this_endPoint = fRoughEndPoint;
1342  }
1343  if(verbose) {
1344  std::cout << "Angle: Start point: (" << this_startPoint.w
1345  << ", " << this_startPoint.t << ")\n";
1346  std::cout << "Angle: End point : (" << this_endPoint.w
1347  << ", " << this_endPoint.t << ")\n";
1348  }
1349  double endStartDiff_x = (this_endPoint.w - this_startPoint.w);
1350  double endStartDiff_y = (this_endPoint.t - this_startPoint.t);
1351  double rms_forward = 0;
1352  double rms_backward = 0;
1353  double mean_forward = 0;
1354  double mean_backward = 0;
1355  //double weight_total = 0;
1356  double hit_counter_forward = 0;
1357  double hit_counter_backward = 0;
1358 
1359  if (verbose && endStartDiff_y == 0 && endStartDiff_x == 0) {
1360  std::cerr << "Error: end point and start point are the same!\n";
1361 
1362  fTimeRecord_ProcName.push_back("RefineDirection");
1363  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1364  return;
1365  }
1366 
1367  double percentage = 0.90;
1368  double percentage_HC = 0.90*fParams.N_Hits_HC/fParams.N_Hits;
1369  const int NBINS=200;
1370  const double wgt = 1.0/fParams.N_Hits;
1371 
1372  // Containers for the angle histograms
1373  std::vector<float> opening_angle_bin(NBINS,0.0 ) ;
1374  std::vector<float> closing_angle_bin(NBINS,0.0 ) ;
1375  std::vector<float> opening_angle_highcharge_bin(NBINS,0.0 ) ;
1376  std::vector<float> closing_angle_highcharge_bin(NBINS,0.0 ) ;
1377  std::vector<float> opening_angle_chargeWgt_bin(NBINS,0.0 ) ;
1378  std::vector<float> closing_angle_chargeWgt_bin(NBINS,0.0 ) ;
1379  //hard coding this for now, should use SetRefineDirectionQMin function
1380  fQMinRefDir = 25;
1381 
1382  int index = -1;
1383  for(auto& hit : fHitVector) {
1384  index ++;
1385  //skip this hit if below minimum cutoff param
1386  if(hit.charge < fQMinRefDir) continue;
1387 
1388  //weight_total = hit.charge;
1389 
1390  // Compute forward mean
1391  double hitStartDiff_x = (hit.w - this_startPoint.w) ;
1392  double hitStartDiff_y = (hit.t - this_startPoint.t) ;
1393 
1394  if (hitStartDiff_x == 0 && hitStartDiff_y == 0) continue;
1395 
1396  double cosangle_start = (endStartDiff_x*hitStartDiff_x +
1397  endStartDiff_y*hitStartDiff_y);
1398  cosangle_start /= ( pow(pow(endStartDiff_x,2)+pow(endStartDiff_y,2),0.5)
1399  * pow(pow(hitStartDiff_x,2)+pow(hitStartDiff_y,2),0.5));
1400 
1401  if(cosangle_start>0) {
1402  // Only take into account for hits that are in "front"
1403  //no weighted average, works better as flat average w/ min charge cut
1404  mean_forward += cosangle_start;
1405  rms_forward += pow(cosangle_start,2);
1406  hit_counter_forward++;
1407  }
1408 
1409  // Compute backward mean
1410  double hitEndDiff_x = (hit.w - this_endPoint.w);
1411  double hitEndDiff_y = (hit.t - this_endPoint.t);
1412  if (hitEndDiff_x == 0 && hitEndDiff_y == 0) continue;
1413 
1414  double cosangle_end = (endStartDiff_x*hitEndDiff_x +
1415  endStartDiff_y*hitEndDiff_y) * (-1.);
1416  cosangle_end /= ( pow(pow(endStartDiff_x,2)+pow(endStartDiff_y,2),0.5)
1417  * pow(pow(hitEndDiff_x,2)+pow(hitEndDiff_y,2),0.5));
1418 
1419  if(cosangle_end>0) {
1420  //no weighted average, works better as flat average w/ min charge cut
1421  mean_backward += cosangle_end;
1422  rms_backward += pow(cosangle_end,2);
1423  hit_counter_backward++;
1424  }
1425 
1426  int N_bins_OPEN = (NBINS-1) * acos(cosangle_start)/PI;
1427  int N_bins_CLOSE = (NBINS-1) * acos(cosangle_end)/PI;
1428  if (N_bins_OPEN < 0) N_bins_OPEN = 0;
1429  if (N_bins_CLOSE < 0) N_bins_CLOSE = 0;
1430  // std::cout << "endStartDiff_x :" << endStartDiff_x << std::endl;
1431  // std::cout << "endStartDiff_y :" << endStartDiff_y << std::endl;
1432  // std::cout << "cosangle_start :" << cosangle_start << std::endl;
1433  // std::cout << "cosangle_end :" << cosangle_end << std::endl;
1434  // std::cout << "N_bins_OPEN :" << N_bins_OPEN << std::endl;
1435  // std::cout << "N_bins_CLOSE :" << N_bins_CLOSE << std::endl;
1436 
1437  opening_angle_chargeWgt_bin[N_bins_OPEN ]
1438  += hit.charge/fParams.sum_charge;
1439  closing_angle_chargeWgt_bin[N_bins_CLOSE]
1440  += hit.charge/fParams.sum_charge;
1441  opening_angle_bin[N_bins_OPEN] += wgt;
1442  closing_angle_bin[N_bins_CLOSE] += wgt;
1443 
1444  //Also fill bins for high charge hits
1445  if(hit.charge > fParams.mean_charge){
1446  opening_angle_highcharge_bin[N_bins_OPEN] += wgt;
1447  closing_angle_highcharge_bin[N_bins_CLOSE] += wgt;
1448  }
1449 
1450  }
1451 
1452  //initialize iterators for the angles
1453  int iBin(0), jBin(0), kBin(0), lBin(0), mBin(0), nBin(0);
1454  double percentage_OPEN(0.0),
1455  percentage_CLOSE(0.0),
1456  percentage_OPEN_HC(0.0),
1457  percentage_CLOSE_HC(0.0), //The last 2 are for High Charge (HC)
1458  percentage_OPEN_CHARGEWGT(0.0),
1459  percentage_CLOSE_CHARGEWGT(0.0);
1460 
1461  for(iBin = 0; percentage_OPEN<= percentage && iBin < NBINS; iBin++)
1462  {
1463  percentage_OPEN += opening_angle_bin[iBin];
1464  }
1465 
1466  for(jBin = 0; percentage_CLOSE<= percentage && jBin < NBINS; jBin++)
1467  {
1468  percentage_CLOSE += closing_angle_bin[jBin];
1469  }
1470 
1471  for(kBin = 0; percentage_OPEN_HC<= percentage_HC && kBin < NBINS; kBin++)
1472  {
1473  percentage_OPEN_HC += opening_angle_highcharge_bin[kBin];
1474  }
1475 
1476  for(lBin = 0; percentage_CLOSE_HC<= percentage_HC && lBin < NBINS; lBin++)
1477  {
1478  percentage_CLOSE_HC += closing_angle_highcharge_bin[lBin];
1479  }
1480 
1481  for(mBin = 0; percentage_OPEN_CHARGEWGT<= percentage && mBin < NBINS; mBin++)
1482  {
1483  percentage_OPEN_CHARGEWGT += opening_angle_chargeWgt_bin[mBin];
1484  }
1485 
1486  for(nBin = 0; percentage_CLOSE_CHARGEWGT<= percentage && nBin < NBINS; nBin++)
1487  {
1488  percentage_CLOSE_CHARGEWGT += closing_angle_chargeWgt_bin[nBin];
1489  }
1490 
1491  double opening_angle = iBin * PI /NBINS ;
1492  double closing_angle = jBin * PI /NBINS ;
1493  double opening_angle_highcharge = kBin * PI /NBINS ;
1494  double closing_angle_highcharge = lBin * PI /NBINS ;
1495  double opening_angle_charge_wgt = mBin * PI /NBINS ;
1496  double closing_angle_charge_wgt = nBin * PI /NBINS ;
1497 
1498  // std::cout<<"opening angle: "<<opening_angle<<std::endl;
1499  // std::cout<<"closing angle: "<<closing_angle<<std::endl;
1500  // std::cout<<"opening high charge angle: "<<opening_angle_highcharge<<std::endl;
1501  // std::cout<<"closing high charge angle: "<<closing_angle_highcharge<<std::endl;
1502  // std::cout<<"opening high charge wgt : "<<opening_angle_charge_wgt<<std::endl;
1503  // std::cout<<"closing high charge wgt : "<<closing_angle_charge_wgt<<std::endl;
1504  // std::cout<<"fCoarseChargeProfile : "<<fCoarseChargeProfile[0]
1505  // << ", " << fCoarseChargeProfile[1] << std::endl;
1506 
1507  double value_1 = closing_angle/opening_angle -1;
1508  // if (fCoarseChargeProfile.size() != 2) std::cout <<" !!!!!!!!!!!!!!!!!! \n\n\n\n";
1509  double value_2 = (fCoarseChargeProfile[0]/fCoarseChargeProfile[1]);
1510  if (value_2 < 100.0 && value_2 > 0.01)
1511  value_2 = log(value_2);
1512  else
1513  value_2 = 0.0;
1514  double value_3 = closing_angle_charge_wgt/opening_angle_charge_wgt -1;
1515 
1516  // if (value_1 > 1.0) value_1 = -1.0/value_1;
1517  // if (value_2 < 1.0) value_2 = -1.0/value_2;
1518  // if (value_3 > 1.0) value_3 = -1.0/value_3;
1519 
1520  // std::cout << "value_1 : " << value_1 << std::endl;
1521  // std::cout << "value_2 : " << value_2 << std::endl;
1522  // std::cout << "value_3 : " << value_3 << std::endl;
1523 
1524  // Using a sigmoid function to determine flipping.
1525  // I am going to weigh all of the values above (1, 2, 3) equally.
1526  // But, since value 2 in particular, I'm going to cut things off at 5.
1527 
1528  if (value_1 > 3 ) value_1 = 3.0;
1529  if (value_1 < -3 ) value_1 = -3.0;
1530  if (value_2 > 3 ) value_2 = 3.0;
1531  if (value_2 < -3 ) value_2 = -3.0;
1532  if (value_3 > 3 ) value_3 = 3.0;
1533  if (value_3 < -3 ) value_3 = -3.0;
1534 
1535  double Exp = exp(value_1 + value_2 + value_3);
1536  double sigmoid = (Exp - 1)/(Exp + 1);
1537 
1538  // std::cout << "sigmoid: " << sigmoid << std::endl;
1539 
1540  bool flip = false;
1541  if (sigmoid < 0) flip = true;
1542  if (flip){
1543  if(verbose) std::cout << "Flipping!" << std::endl;
1544  std::swap(opening_angle,closing_angle);
1545  std::swap(opening_angle_highcharge,closing_angle_highcharge);
1546  std::swap(opening_angle_charge_wgt,closing_angle_charge_wgt);
1549  }
1550  else if(verbose){
1551  std::cout << "Not Flipping!\n";
1552  }
1553 
1554  fParams.opening_angle = opening_angle;
1555  fParams.opening_angle_charge_wgt = opening_angle_charge_wgt;
1556  fParams.closing_angle = closing_angle;
1557  fParams.closing_angle_charge_wgt = closing_angle_charge_wgt;
1558 
1559  fFinishedRefineDirection = true;
1560 
1561  fTimeRecord_ProcName.push_back("RefineDirection");
1562  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1563 
1564  // return;
1565  } //end RefineDirection
double closing_angle_charge_wgt
Same for charge_wgt.
Definition: ClusterParams.h:44
std::vector< std::string > fTimeRecord_ProcName
std::vector< util::PxHit > fHitVector
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
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
std::vector< double > fCoarseChargeProfile
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
#define PI
double closing_angle
Width of angular distubtion wrt endpoint.
Definition: ClusterParams.h:43
void cluster::ClusterParamsAlg::RefineStartPointAndDirection ( bool  override = false)

Definition at line 1596 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::direction, cluster::cluster_params::end_point, fFinishedGetProfileInfo, fFinishedRefineStartPointAndDirection, fParams, fRoughBeginPoint, fRoughEndPoint, fTimeRecord_ProcName, fTimeRecord_ProcTime, GetProfileInfo(), RefineDirection(), RefineStartPoints(), cluster::cluster_params::start_point, util::PxPoint::t, verbose, and util::PxPoint::w.

Referenced by FillParams(), and GetParams().

1596  {
1597  // This function is meant to pick the direction.
1598  // It refines both the start and end point, and then asks
1599  // if it should flip.
1600 
1601  TStopwatch localWatch;
1602  localWatch.Start();
1603 
1604  if(verbose) std::cout << " here!!! " << std::endl;
1605 
1606  if(!override) { //Override being set, we skip all this logic.
1607  //OK, no override. Stop if we're already finshed.
1609  fTimeRecord_ProcName.push_back("RefineStartPointAndDirection");
1610  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1611  return;
1612  }
1613  //Try to run the previous function if not yet done.
1615  } else {
1616  //Try to run the previous function if not yet done.
1618  }
1619  if(verbose){
1620  std::cout << "REFINING .... " << std::endl;
1621  std::cout << " Rough start and end point: " << std::endl;
1622  std::cout << " s: (" << fParams.start_point.w << ", "
1623  << fParams.start_point.t << ")" << std::endl;
1624  std::cout << " e: (" << fParams.end_point.w << ", "
1625  << fParams.end_point.t << ")" << std::endl;
1626  }
1628  if(verbose){
1629  std::cout << " Once Refined start and end point: " << std::endl;
1630  std::cout << " s: (" << fParams.start_point.w << ", "
1631  << fParams.start_point.t << ")" << std::endl;
1632  std::cout << " e: (" << fParams.end_point.w << ", "
1633  << fParams.end_point.t << ")" << std::endl;
1636  }
1638  if(verbose) {
1639  std::cout << " Twice Refined start and end point: " << std::endl;
1640  std::cout << " s: (" << fParams.start_point.w << ", "
1641  << fParams.start_point.t << ")" << std::endl;
1642  std::cout << " e: (" << fParams.end_point.w << ", "
1643  << fParams.end_point.t << ")" << std::endl;
1646  }
1647  RefineDirection();
1648  if(verbose) {
1649  std::cout << " Final start and end point: " << std::endl;
1650  std::cout << " s: (" << fParams.start_point.w << ", "
1651  << fParams.start_point.t << ")" << std::endl;
1652  std::cout << " e: (" << fParams.end_point.w << ", "
1653  << fParams.end_point.t << ")" << std::endl;
1654  }
1656 
1657  fTimeRecord_ProcName.push_back("RefineStartPointAndDirection");
1658  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1659  return;
1660  }
void GetProfileInfo(bool override=false)
std::vector< std::string > fTimeRecord_ProcName
void RefineStartPoints(bool override=false)
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
double w
Definition: PxUtils.h:10
void RefineDirection(bool override=false)
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
util::PxPoint end_point
end point
Definition: ClusterParams.h:25
void cluster::ClusterParamsAlg::RefineStartPoints ( bool  override = false)

Calculates the following variables: length width

Parameters
override[description]

Calculates the following variables: length width hit_density_1D hit_density_2D direction

Parameters
override[description]

Definition at line 869 of file ClusterParamsAlg.cxx.

References util::PxHit::charge, d, cluster::cluster_params::end_point, fFinishedRefineStartPoints, fGSer, fHitVector, fParams, fRough2DSlope, fRoughBeginPoint, fRoughEndPoint, fTimeRecord_ProcName, fTimeRecord_ProcTime, util::GeometryUtilities::Get2DDistance(), cluster::cluster_params::hit_density_1D, cluster::cluster_params::hit_density_2D, cluster::cluster_params::length, n, cluster::cluster_params::N_Hits, util::PxPoint::plane, util::GeometryUtilities::SelectLocalHitlist(), cluster::cluster_params::start_point, util::PxPoint::t, verbose, util::PxPoint::w, cluster::cluster_params::width, and z.

Referenced by GetFinalSlope(), GetParams(), MultipleHitDensity(), and RefineStartPointAndDirection().

869  {
870 
871  TStopwatch localWatch;
872  localWatch.Start();
873 
874  //
875  // Why override is not used here? Kazu 05/01/2014
876  // Put a stupid line here to avoid compiler warning
877  if(!override) override = true;
878 
879  // if(!override) { //Override being set, we skip all this logic.
880  // //OK, no override. Stop if we're already finshed.
881  // if (fFinishedRefineStartPoints) return;
882  // //Try to run the previous function if not yet done.
883  // if (!fFinishedRefineDirection) RefineDirection(true);
884  // } else {
885  // if (!fFinishedRefineDirection) RefineDirection(true);
886  // }
887 
888 
889  // need to define physical direction with openind angles and pass that to Ryan's line finder.
890 
891 
892  // Direction decision has been moved entirely to Refine Direction()
893  // opening and closing angles are already set
894  // they are associated with the start and endpoints.
895  // If refine direction opted to switch the shower direction,
896  // it also switched opening and closing angles.
897 
899  // fParams.direction= ...
900 
901  // Direction is decided by RefineDirection()
902  util::PxPoint startHit,endHit;
903  startHit=fRoughBeginPoint;
904  endHit=fRoughEndPoint;
905 
906 
908  //Ryan's Shower Strip finder work here.
909  //First we need to define the strip width that we want
910  double d=0.6;//this is the width of the strip.... this needs to be tuned to something.
911  //===============================================================================================================
912  // Will need to feed in the set of hits that we want.
913  // const std::vector<util::PxHit*> whole;
914  std::vector <const util::PxHit*> subhit;
915  double linearlimit=8;
916  double ortlimit=12;
917  double lineslopetest(fRough2DSlope);
918  util::PxHit averageHit;
919  //also are we sure this is actually doing what it is supposed to???
920  //are we sure this works?
921  //is anybody even listening? Were they eaten by a grue?
922  fGSer->SelectLocalHitlist(fHitVector,subhit, startHit,
923  linearlimit,ortlimit,lineslopetest,
924  averageHit);
925 
926 
927  if(!(subhit.size()) || subhit.size()<3) {
928  if(verbose) std::cout<<"Subhit list is empty or too small. Using rough start/end points..."<<std::endl;
929  // GetOpeningAngle();
932  // fRoughEndPoint
933  // fRoughEndPoint
934  // and use them to get the axis
935 
937 
938  fTimeRecord_ProcName.push_back("RefineStartPoints");
939  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
940 
941  return;
942  }
943 
944  double avgwire= averageHit.w;
945  double avgtime= averageHit.t;
946  //vertex in tilda-space pair(x-til,y-til)
947  std::vector<std::pair<double,double>> vertil;
948  // vertil.clear();// This isn't needed?
949  vertil.reserve(subhit.size() * subhit.size());
950  //vector of the sum of the distance of a vector to every vertex in tilda-space
951  std::vector<double> vs;
952  // vs.clear();// this isn't needed?
953  // $$This needs to be corrected//this is the good hits that are between strip
954  std::vector<const util::PxHit*> ghits;
955  ghits.reserve(subhit.size());
956  int n=0;
957  double fardistcurrent=0;
958  util::PxHit startpoint;
959  double gwiretime=0;
960  double gwire=0;
961  double gtime=0;
962  double gwirewire=0;
963  //KAZU!!! I NEED this too//this will need to come from somewhere...
964  //This is some hit that is from the way far side of the entire shower cluster...
966 
967  //=============building the local vector===================
968  //this is clumsy... but just want to get something to work for now.
969  //THIS is REAL Horrible and CLUMSY... We can make things into helper functions soon.
970  std::vector<util::PxHit> returnhits;
971  std::vector<double> radiusofhit;
972  std::vector<int> closehits;
973  //unsigned int minhits=0;
974  //double maxreset=0;
975  // double avgwire=0;
976  // double avgtime=0;
977  // if(minhits<fHitVector.size()){
978  // for(auto & hit : fHitVector){
979  // std::pair<double,double> rv(fRoughEndPoint.w,fRoughEndPoint.t);
980  // closehits.clear();
981  // radiusofhit.clear();
982  // returnhits.clear();
983  // for( unsigned int a=0; a<hit.size(); a++){
984  // double d= sqrt( pow((rv.first-hit.w),2) + pow((rv.second-hit.t),2) );
985  // maxreset+=d;
986  // radiusofhit.push_back(d);}
987  // for(unsigned int b=0; b<minhits; b++){
988  // int minss= std::min_element(radiusofhit.begin(),radiusofhit.end())-radiusofhit.begin();
989  // closehits.push_back(minss);
990  // radiusofhit[minss] = maxreset;}
991  //now make the vector o just the close hits using the closehit index
992  // for( unsigned int k=0; k < closehits.size(); k++){
993  // //first find the average wire and time for each set of hits... and make that the new origin before the transpose.
994  // avgwire+=fHitVector[closehits[k]].w;
995  // avgtime+=fHitVector[closehits[k]].t;
996  // returnhits.push_back(fHitVector[closehits[k]]);}
997  // }//if fHitVector is big enough
998 
999  // avgwire=avgwire/closehits.size();
1000  // avgtime=avgtime/closehits.size();
1001  // subhit=returnhits;
1002 
1003  //==============================================================================
1004  //Now we need to do the transformation into "tilda-space"
1005  for(unsigned int a=0; a<subhit.size();a++){
1006  for(unsigned int b=a+1;b<subhit.size();b++){
1007  if(subhit.at(a)->w != subhit.at(b)->w){
1008  double xtil = ((subhit.at(a)->t - avgtime) - (subhit.at(b)->t - avgtime));
1009  xtil /= ((subhit.at(a)->w - avgwire)-(subhit.at(b)->w - avgwire));
1010  double ytil = (subhit.at(a)->w - avgwire)*xtil -(subhit.at(a)->t - avgtime);
1011  //push back the tilda vertex point on the pair
1012  std::pair<double,double> tv(xtil,ytil);
1013  vertil.push_back(tv);
1014  }//if Wires are not the same
1015  }//for over b
1016  }//for over a
1017  // look at the distance from a tilda-vertex to all other tilda-verticies
1018  for(unsigned int z=0;z<vertil.size();z++){
1019  double vd=0;//the distance for vertex... just needs to be something 0
1020  for(unsigned int c=0; c<vertil.size();c++)
1021 
1022  vd+= sqrt(pow((vertil.at(z).first - vertil.at(c).first),2)
1023  + pow((vertil.at(z).second-vertil.at(c).second),2));
1024  vs.push_back(vd);
1025  vd=0.0;
1026  }//for z loop
1027 
1028  if(vs.size()==0) //al hits on same wire?!
1029  {
1030  if(verbose) std::cout<<"vertil list is empty. all subhits are on the same wire?"<<std::endl;
1031  // GetOpeningAngle();
1034  // fRoughEndPoint
1035  // fRoughEndPoint
1036  // and use them to get the axis
1037 
1039 
1040  fTimeRecord_ProcName.push_back("RefineStartPoints");
1041  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1042  return;
1043  }
1044  //need to find the min of the distance of vertex in tilda-space
1045  //this will get me the area where things are most linear
1046  int minvs= std::min_element(vs.begin(),vs.end())-vs.begin();
1047  // now use the min position to find the vertex in tilda-space
1048  //now need to look a which hits are between the tilda lines from the points
1049  //in the tilda space everything in wire time is based on the new origin which is at the average wire/time
1050  double tilwire=vertil.at(minvs).first;//slope
1051  double tiltimet=-vertil.at(minvs).second+d*sqrt(1+pow(tilwire,2));//negative cept is accounted for top strip
1052  double tiltimeb=-vertil.at(minvs).second-d*sqrt(1+pow(tilwire,2));//negative cept is accounted for bottom strip
1053  // look over the subhit list and ask for which are inside of the strip
1054  for(unsigned int a=0; a<subhit.size(); a++){
1055  double dtstrip= (-tilwire * (subhit.at(a)->w - avgwire)
1056  +(subhit.at(a)->t - avgtime)-tiltimet)/sqrt(tilwire*tilwire+1);
1057  double dbstrip= (-tilwire * (subhit.at(a)->w - avgwire)
1058  +(subhit.at(a)->t - avgtime)-tiltimeb)/sqrt(tilwire*tilwire+1);
1059 
1060  if((dtstrip<0.0 && dbstrip>0.0)||(dbstrip<0.0 && dtstrip>0.0)){
1061  ghits.push_back(subhit.at(a));
1062  }//if between strip
1063  }//for a loop
1064  //=======================Do stuff with the good hits============================
1065 
1066  //of these small set of hits just fit a simple line.
1067  //then we will need to see which of these hits is farthest away
1068 
1069  for(unsigned int g=0; g<ghits.size();g++){
1070  // should call the helper funtion to do the fit
1071  //but for now since there is no helper function I will just write it here......again
1072  n+=1;
1073  gwiretime+= ghits.at(g)->w * ghits.at(g)->t;
1074  gwire+= ghits.at(g)->w;
1075  gtime+= ghits.at(g)->t;
1076  gwirewire+= ghits.at(g)->w * ghits.at(g)->w;
1077  // now work on calculating the distance in wire time space from the far point
1078  //farhit needs to be a hit that is given to me
1079  double fardist= sqrt(pow(ghits.at(g)->w - farhit.w,2)+pow(ghits.at(g)->t - farhit.t,2));
1080  //come back to this... there is a better way to do this probably in the loop
1081  //there should also be a check that the hit that is farthest away has subsequent hits after it on a few wires
1082  if(fardist>fardistcurrent){
1083  fardistcurrent=fardist;
1084  //if fardist... this is the point to use for the start point
1085  startpoint.t = ghits.at(g)->t;
1086  startpoint.w = ghits.at(g)->w;
1087  startpoint.plane = ghits.at(g)->plane;
1088  startpoint.charge = ghits.at(g)->charge;
1089  }
1090  }//for ghits loop
1091  //This can be the new start point
1092  // std::cout<<"Here Kazu"<<std::endl;
1093  // std::cout<<"Ryan This is the new SP ("<<startpoint.w<<" , "<<startpoint.t<<")"<<std::endl;
1094  // double gslope=(n* gwiretime- gwire*gtime)/(n*gwirewire -gwire*gwire);
1095  // double gcept= gtime/n -gslope*(gwire/n);
1096 
1097  //should this be here? Id argue this might be done once outside.
1098  fParams.length=fGSer->Get2DDistance((util::PxPoint *)&startpoint,&endHit);
1099  if(fParams.length)
1101  else
1103 
1104  if(fParams.length && fParams.width)
1106  else
1108 
1109 
1110  fParams.start_point=startpoint;
1111 
1112  static int count(0);
1113  count ++;
1114  // std::cout << "Completed refine start point " << count << " times!\n";
1115 
1117  // Report();
1118 
1119  fTimeRecord_ProcName.push_back("RefineStartPoints");
1120  fTimeRecord_ProcTime.push_back(localWatch.RealTime());
1121 
1122  return;
1123  }
util::GeometryUtilities * fGSer
std::vector< std::string > fTimeRecord_ProcName
Double_t Get2DDistance(Double_t wire1, Double_t time1, Double_t wire2, Double_t time2) const
Double_t z
Definition: plot.C:279
std::vector< util::PxHit > fHitVector
cluster::cluster_params fParams
double t
Definition: PxUtils.h:11
Float_t d
Definition: plot.C:237
std::vector< double > fTimeRecord_ProcTime
util::PxPoint start_point
start point
Definition: ClusterParams.h:24
void SelectLocalHitlist(const std::vector< util::PxHit > &hitlist, std::vector< const util::PxHit * > &hitlistlocal, util::PxPoint &startHit, Double_t &linearlimit, Double_t &ortlimit, Double_t &lineslopetest)
double w
Definition: PxUtils.h:10
util::PxPoint end_point
end point
Definition: ClusterParams.h:25
double charge
area charge
Definition: PxUtils.h:39
Char_t n[5]
unsigned int plane
Definition: PxUtils.h:12
template<typename Stream >
void cluster::ClusterParamsAlg::Report ( Stream &  stream) const

Definition at line 243 of file ClusterParamsAlg.cxx.

References fFinishedGetAverages, fFinishedGetEndCharges, fFinishedGetFinalSlope, fFinishedGetProfileInfo, fFinishedGetRoughAxis, fFinishedRefineDirection, fFinishedRefineStartPoints, fParams, cluster::cluster_params::Report(), and verbose.

Referenced by SetVerbose().

243  {
244  if(verbose) {
245  stream << "ClusterParamsAlg Report: " << "\n"
246  << "\tFinishedGetAverages " << fFinishedGetAverages << "\n"
247  << "\tFinishedGetRoughAxis " << fFinishedGetRoughAxis << "\n"
248  << "\tFinishedGetProfileInfo " << fFinishedGetProfileInfo << "\n"
249  << "\tFinishedRefineStartPoints " << fFinishedRefineStartPoints << "\n"
250  << "\tFinishedRefineDirection " << fFinishedRefineDirection << "\n"
251  << "\tFinishedGetFinalSlope " << fFinishedGetFinalSlope << "\n"
252  << "\tFinishedGetEndCharges " << fFinishedGetEndCharges << "\n"
253  << "--------------------------------------" << "\n";
254  fParams.Report(stream);
255  }
256  }
void Report(Stream &os) const
cluster::cluster_params fParams
const util::PxPoint& cluster::ClusterParamsAlg::RoughEndPoint ( )
inline

Definition at line 202 of file ClusterParamsAlg.h.

References fRoughEndPoint.

202 {return fRoughEndPoint;}
double cluster::ClusterParamsAlg::RoughIntercept ( )
inline
double cluster::ClusterParamsAlg::RoughSlope ( )
inline

Definition at line 204 of file ClusterParamsAlg.h.

References fRough2DSlope.

const util::PxPoint& cluster::ClusterParamsAlg::RoughStartPoint ( )
inline

Definition at line 201 of file ClusterParamsAlg.h.

References fRoughBeginPoint.

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

Definition at line 52 of file ClusterParamsAlg.cxx.

References fHitVector, fMinNHits, fPlane, Initialize(), and verbose.

Referenced by ClusterParamsAlg(), MinNHits(), and cluster::StandardClusterParamsAlg::SetHits().

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

Definition at line 61 of file ClusterParamsAlg.h.

References fMinNHits.

Referenced by cmtool::CMManagerBase::SetClusters().

61 { fMinNHits = nhit; }
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 195 of file ClusterParamsAlg.h.

References FillPolygon(), fNeuralNetPath, GetOpeningAngle(), and s.

195 {fNeuralNetPath = s;}
Float_t s
Definition: plot.C:23
void cluster::ClusterParamsAlg::SetPlane ( int  p)

Definition at line 118 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::end_point, fHitVector, fParams, fPlane, fRoughBeginPoint, fRoughEndPoint, util::PxPoint::plane, and cluster::cluster_params::start_point.

Referenced by Plane().

118  {
119  fPlane = p;
120  for(auto& h : fHitVector) h.plane = p;
123  }
std::vector< util::PxHit > fHitVector
cluster::cluster_params fParams
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 69 of file ClusterParamsAlg.h.

References fQMinRefDir.

69 { fQMinRefDir = qmin; }
void cluster::ClusterParamsAlg::SetVerbose ( bool  yes = true)
inline
double cluster::ClusterParamsAlg::StartCharge ( 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 1770 of file ClusterParamsAlg.cxx.

References EndCharge(), fHitVector, and IntegrateFitCharge().

Referenced by GetEndCharges(), IntegrateFitCharge(), RoughIntercept(), and cluster::StandardClusterParamsAlg::StartCharge().

1771  {
1772  switch (fHitVector.size()) {
1773  case 0: return 0.;
1774  case 1: return fHitVector.back().charge;
1775  // the "default" is the rest of the function
1776  } // switch
1777 
1778  // this is the range of the fit:
1779  const unsigned int fit_first_bin = 0, fit_last_bin = nbins;
1780 
1781  return IntegrateFitCharge(0., length, fit_first_bin, fit_last_bin);
1782  } // ClusterParamsAlg::StartCharge()
std::vector< util::PxHit > fHitVector
double IntegrateFitCharge(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.
template<typename Stream >
void cluster::ClusterParamsAlg::TimeReport ( Stream &  stream) const

Definition at line 397 of file ClusterParamsAlg.h.

References fTimeRecord_ProcName, and fTimeRecord_ProcTime.

Referenced by SetVerbose().

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

Definition at line 1662 of file ClusterParamsAlg.cxx.

References enableFANN, GetFANNVector(), and verbose.

Referenced by FillParams(), and GetParams().

1662  {
1663  if(!override) return;
1664  if(verbose) std::cout << " ---- Trying T/S sep. ------ \n";
1665  if (enableFANN){
1666  if(verbose) std::cout << " ---- Doing T/S sep. ------- \n";
1667  std::vector<float> FeatureVector, outputVector;
1668  GetFANNVector(FeatureVector);
1669 // ::cluster::FANNService::GetME()->GetFANNModule().run(FeatureVector,outputVector);
1670  // fParams.trackness = outputVector[1];
1671  // fParams.showerness = outputVector[0];
1672  }
1673  else{
1674  if(verbose) std::cout << " ---- Failed T/S sep. ------ \n";
1675  }
1676  }
void GetFANNVector(std::vector< float > &data)

Member Data Documentation

bool cluster::ClusterParamsAlg::enableFANN
protected
double cluster::ClusterParamsAlg::fBeginIntercept
protected

Definition at line 316 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo().

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

Definition at line 293 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), and Initialize().

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

Definition at line 299 of file ClusterParamsAlg.h.

Referenced by EndCharge(), GetFinalSlope(), GetProfileInfo(), and IntegrateFitCharge().

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

Definition at line 302 of file ClusterParamsAlg.h.

Referenced by GetFinalSlope().

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

Definition at line 300 of file ClusterParamsAlg.h.

Referenced by GetFinalSlope(), GetProfileInfo(), and RefineDirection().

int cluster::ClusterParamsAlg::fCoarseNbins
protected

Definition at line 306 of file ClusterParamsAlg.h.

Referenced by GetFinalSlope(), and GetProfileInfo().

double cluster::ClusterParamsAlg::fEndIntercept
protected

Definition at line 317 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo().

bool cluster::ClusterParamsAlg::fFinishedGetAverages
protected

Definition at line 322 of file ClusterParamsAlg.h.

Referenced by GetAverages(), GetRoughAxis(), Initialize(), and Report().

bool cluster::ClusterParamsAlg::fFinishedGetEndCharges
protected

Definition at line 330 of file ClusterParamsAlg.h.

Referenced by GetEndCharges(), Initialize(), and Report().

bool cluster::ClusterParamsAlg::fFinishedGetFinalSlope
protected

Definition at line 327 of file ClusterParamsAlg.h.

Referenced by GetFinalSlope(), Initialize(), and Report().

bool cluster::ClusterParamsAlg::fFinishedGetProfileInfo
protected
bool cluster::ClusterParamsAlg::fFinishedGetRoughAxis
protected

Definition at line 323 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), GetRoughAxis(), Initialize(), and Report().

bool cluster::ClusterParamsAlg::fFinishedRefineDirection
protected

Definition at line 326 of file ClusterParamsAlg.h.

Referenced by Initialize(), RefineDirection(), and Report().

bool cluster::ClusterParamsAlg::fFinishedRefineStartPointAndDirection
protected

Definition at line 328 of file ClusterParamsAlg.h.

Referenced by Initialize(), and RefineStartPointAndDirection().

bool cluster::ClusterParamsAlg::fFinishedRefineStartPoints
protected
bool cluster::ClusterParamsAlg::fFinishedTrackShowerSep
protected

Definition at line 329 of file ClusterParamsAlg.h.

Referenced by Initialize().

util::GeometryUtilities* cluster::ClusterParamsAlg::fGSer
protected
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 286 of file ClusterParamsAlg.h.

Referenced by EndCharge(), FillPolygon(), GetAverages(), GetFinalSlope(), GetHitVector(), GetNHits(), GetProfileInfo(), GetRoughAxis(), Initialize(), MultipleHitDensity(), MultipleHitWires(), RefineDirection(), RefineStartPoints(), SetHits(), SetPlane(), and StartCharge().

double cluster::ClusterParamsAlg::fInterHigh_side
protected

Definition at line 318 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo().

double cluster::ClusterParamsAlg::fInterLow_side
protected

Definition at line 319 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo().

size_t cluster::ClusterParamsAlg::fMinNHits
protected

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

Definition at line 280 of file ClusterParamsAlg.h.

Referenced by ClusterParamsAlg(), MinNHits(), SetHits(), and SetMinNHits().

std::string cluster::ClusterParamsAlg::fNeuralNetPath

Definition at line 378 of file ClusterParamsAlg.h.

Referenced by setNeuralNetPath().

int cluster::ClusterParamsAlg::fPlane
protected

Definition at line 294 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), Plane(), SetHits(), and SetPlane().

double cluster::ClusterParamsAlg::fProfileIntegralBackward
protected

Definition at line 310 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), and Initialize().

double cluster::ClusterParamsAlg::fProfileIntegralForward
protected

Definition at line 309 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), and Initialize().

int cluster::ClusterParamsAlg::fProfileMaximumBin
protected

Definition at line 308 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), and Initialize().

int cluster::ClusterParamsAlg::fProfileNbins
protected

Definition at line 307 of file ClusterParamsAlg.h.

Referenced by GetFinalSlope(), and GetProfileInfo().

double cluster::ClusterParamsAlg::fProjectedLength
protected

Definition at line 311 of file ClusterParamsAlg.h.

Referenced by EndCharge(), GetFinalSlope(), GetProfileInfo(), and IntegrateFitCharge().

double cluster::ClusterParamsAlg::fQMinRefDir
protected

Definition at line 297 of file ClusterParamsAlg.h.

Referenced by RefineDirection(), and SetRefineDirectionQMin().

double cluster::ClusterParamsAlg::fRough2DIntercept
protected
double cluster::ClusterParamsAlg::fRough2DSlope
protected
util::PxPoint cluster::ClusterParamsAlg::fRoughBeginPoint
protected
util::PxPoint cluster::ClusterParamsAlg::fRoughEndPoint
protected
std::vector<std::string> cluster::ClusterParamsAlg::fTimeRecord_ProcName
std::vector<double> cluster::ClusterParamsAlg::fTimeRecord_ProcTime

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