9 #include "Utilities/func/MathUtil.h" 37 for (
unsigned int n=0;
n<pchits.size(); ++
n){
49 default:
assert(0 &&
"Unknown view");
60 default:
assert(0 &&
"Unknown view");
88 std::vector<art::Ptr<caldp::PCHit> > all;
123 for(
unsigned int i=0;
i<
NPCHit(view); ++
i){
133 for(
unsigned int i=0;
i<
NPCHit(view); ++
i){
143 for(
unsigned int i=0;
i<
NPCHit(view); ++
i){
153 for(
unsigned int i=0;
i<
NPCHit(view); ++
i){
162 double ret = DBL_MAX;
163 for(
unsigned int i=0;
i<
NPCHit(view); ++
i){
172 double ret = -DBL_MAX;
173 for(
unsigned int i=0;
i<
NPCHit(view); ++
i){
182 double ret = DBL_MAX;
183 for(
unsigned int i=0;
i<
NPCHit(view); ++
i){
184 if(
PCHit(view,
i)->FlightLen() <
ret) ret =
PCHit(view,
i)->FlightLen();
192 double ret = -DBL_MAX;
193 for(
unsigned int i=0;
i<
NPCHit(view); ++
i){
194 if(
PCHit(view,
i)->FlightLen() >
ret) ret =
PCHit(view,
i)->FlightLen();
211 double stdev = 144.0;
214 stdev = 231594.0/(2267.16+
pow(pe,2.01011)) + 9.55689;
217 stdev = 848357.0/(5734.75+
pow(pe,2.20060)) + 142.221;
222 stdev = 2229.53/(77.0043+
pow(pe,1.22743)) + 4.89355;
225 stdev = 6136470.0/(185518.0+
pow(pe,2.92360)) + 31.6071;
254 + x*(-3.81279e-23))))))));
279 double fspeed = 15.7446;
281 fspeed = 15.7446/(((15.7446*1.63517)/
d)+1.0);
293 std::vector<double> *resids,
306 std::vector<double>
time;
307 std::vector<double> rDist;
308 std::vector<double> sterror;
309 std::vector<double>
weight;
314 double pe = pch->
PE();
318 tns = pch->
TNS() + tw;
322 sterror.push_back(stdev);
323 weight.push_back(1.0/(stdev*stdev));
325 double slope, intercept = -99.0,
x1,
y1;
326 double fspeed = -99.0;
330 if (rDist.size() < 2)
break;
332 *chisqr =
util::LinFit(rDist, time, weight, slope, intercept);
335 y1 = slope*
x1 + intercept;
337 unsigned int maxpos = 0;
338 for(
unsigned int j=0;
j<rDist.size(); ++
j){
339 double recresid = time[
j] - (fspeed*(rDist[
j]-
x1)+y1);
340 resids->push_back(recresid);
341 if (
fabs(recresid/sterror[
j]) >
fabs(maxresid)){
342 maxresid =
fabs(recresid/sterror[j]);
347 if (rDist.size() == 2)
break;
350 if (
fabs(maxresid) > timeResid){
351 weight.erase(weight.begin()+maxpos);
352 rDist.erase(rDist.begin()+maxpos);
353 time.erase(time.begin()+maxpos);
354 sterror.erase(sterror.begin()+maxpos);
355 resids->erase(resids->begin()+maxpos);
356 if (removeHits ==
true){
375 double fiberVel)
const 385 std::vector<double>
time;
386 std::vector<double> rDist;
387 std::vector<double>
weight;
392 double pe = pch->
PE();
399 tns = pch->
TNS() + tw;
400 time.push_back(tns - pch->
ReadDist()/fiberVel);
404 weight.push_back(1.0/(stdev*stdev));
407 double slope, intercept,
x1,
y1;
408 double mspeed = -5.0;
413 if (rDist.size() < 2)
return -5.0;
414 *chisqr =
util::LinFit(rDist, time, weight, slope, intercept);
417 y1 = slope*x1 + intercept;
418 for(
unsigned int j=0;
j<rDist.size(); ++
j){
419 double recresid = time[
j] - (mspeed*(rDist[
j]-
x1)+y1);
420 resids->push_back(recresid);
double MinTNS(geo::View_t view=geo::kXorY) const
return min pchit time for a specified view
int MaxCell(geo::View_t view) const
return max cell number in the specified view
double CalculateFiberVelocity(double timeResid, std::vector< double > *resids, double *icept, double *chisqr, novadaq::cnv::DetId detectorID, bool removeHits=false, bool overwrite=false)
float TNS() const
Return uncorrected hit time.
float FlightLen() const
Return path from start of track to current cell.
fvar< T > fabs(const fvar< T > &x)
unsigned int NPCHit() const
return total number of hits
Float_t y1[n_points_granero]
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Float_t x1[n_points_granero]
bool GoodTime() const
Return quality of timing fit for cell.
Vertical planes which measure X.
art::Ptr< caldp::PCHit > YPCHit(unsigned int yIdx) const
get the ith y view pchit
::xsd::cxx::tree::time< char, simple_type > time
double MaxTNS(geo::View_t view=geo::kXorY) const
return max pchit time for a specified view
int MaxPlane(geo::View_t view) const
return max plane number in the specified view
Horizontal planes which measure Y.
Far Detector at Ash River, MN.
virtual void Add(const art::Ptr< caldp::PCHit > &pchit)
add a pchit to the track
float PE() const
Return PE value.
std::vector< art::Ptr< caldp::PCHit > > fDropPCHits
contains collection of pchits dropped by fiber fit
void RemovePCHit(const art::Ptr< caldp::PCHit > pch)
remove a pchit from the current cluster
double MinPath(geo::View_t view=geo::kXorY) const
return shortest path from start of track to a hit in a view
Near Detector in the NuMI cavern.
int MinCell(geo::View_t view) const
return min cell number in the specified view
double MaxPath(geo::View_t view=geo::kXorY) const
return longest path from start of track to a hit in a view
geo::View_t View() const
Return view.
unsigned int NXPCHit() const
return number of hits in x view
Histograms used by attenuation calibration.
art::Ptr< caldp::PCHit > PCHit(geo::View_t view, unsigned int viewIdx) const
get the ith pchit from a specified view
art::Ptr< caldp::PCHit > XPCHit(unsigned int xIdx) const
get the ith x view pchit
float ReadDist() const
Return distance to readout, with pigtail.
static const double sCmPerNsec
double CalcFiberVelocity(float d, novadaq::cnv::DetId detectorID) const
Calculate fiber speed as a function of distance from readout.
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
assert(nhit_max >=nhit_nbins)
double ReadoutDistCorrection(float x, novadaq::cnv::DetId detectorID) const
Correct for timewalk effect as a function of distance to readout.
std::vector< art::Ptr< caldp::PCHit > > fYPCHits
double fFiberVelocity
optical fiber velocity for track
double CalculateMuonVelocity(std::vector< double > *resids, double *icept, double *chisqr, novadaq::cnv::DetId detectorID, double fiberVel=9999.0) const
std::vector< art::Ptr< caldp::PCHit > > AllPCHits() const
return the collection of all pchits
int MinPlane(geo::View_t view) const
return min plane number in the specified view
unsigned int NYPCHit() const
return number of hits in y view
std::vector< art::Ptr< caldp::PCHit > > fXPCHits
contain collections of pchits
double TNSUncertainty(float pe, bool goodTime, novadaq::cnv::DetId detectorID) const
Calculate the uncertainty in the time based on pulse height.