37 float ret = trackdir*beamdir;
46 double xyz[3], getx, gety, getz, intercept;
47 std::vector<double> zs, ts, qs;
48 TVector3 tdir = track->
Dir();
49 TVector3 tstart = track->
Start();
50 TVector3 tend = track->
Stop();
55 for(
unsigned int hitIdx = 0; hitIdx < track->
NCell(); ++hitIdx){
69 px = -tdir.
Z()/tdir.X();
75 xyz[0]+50*px, xyz[2]+50*pz, getx, getz);
76 zs.push_back((getz-tstart.Z())/tdir.Z());
82 py = -tdir.Z()/tdir.Y();
88 xyz[1]+50*py, xyz[2]+50*pz, gety, getz);
89 zs.push_back((getz-tstart.Z())/tdir.Z());
92 ts.push_back(rhit.
T());
95 double par1 = -0.377*rhit.
PECorr();
99 double error = 87.3+par1+par2+par3+par4;
100 if(error<1) error = 1;
101 qs.push_back(1/(error*error));
109 slope =
fitslope(zs, ts, qs, chisq, intercept, 1);
122 float CosRejFxs::fitslope(
const std::vector< double > &
l,
const std::vector< double > &
t,
const std::vector< double > &q,
double &chisq,
double &inter,
bool linflag) {
124 assert(t.size()==l.size());
127 double ztStart[2] = {0.};
128 double ztEnd[2] = {0.};
131 &ztStart[0], &ztStart[1],
132 &ztEnd[0], &ztEnd[1]);
134 else if(linflag==
true) {
136 &ztStart[0], &ztStart[1],
137 &ztEnd[0], &ztEnd[1]);
139 const double ztSlope = (ztEnd[1]-ztStart[1])/(ztEnd[0]-ztStart[0]);
141 inter = ztStart[1]-(ztSlope*ztStart[0]);
150 assert(t.size()==l.size());
157 for (
unsigned int i=1;
i<t.size(); ++
i) {
162 double norm = t.size();
167 double tstart = tave - 500;
168 double tend = tave + 500;
171 assert(flag == 0 || flag == 1);
173 lstart = lave - (c * 500);
174 lend = lave + (c * 500);
177 lstart = lave + (c * 500);
178 lend = lave - (c * 500);
181 for (
unsigned int i=0;
i<t.size(); ++
i) {
190 double tscattmax=0, tscatttot=0, tscatt;
196 dir = (traj2-traj1).
Unit();
197 if((dir-olddir).Mag()<1
e-7)
continue;
198 if(
p>0&&dir.Mag()<=1.0&&olddir.Mag()<=1.0) {
199 tscatt =
std::acos(dir.Dot(olddir)) * 180 / TMath::Pi();
200 if(tscatt > tscattmax) tscattmax = tscatt;
213 for (
unsigned int sliceHit = 0; sliceHit < sliceHits.
size(); ++sliceHit){
215 for (
unsigned int trackHit = 0; trackHit < trackHits.
size(); ++trackHit){
216 match = (sliceHits[sliceHit] == trackHits[trackHit]);
219 if (!match) vertexCluster.
Add(sliceHits[sliceHit]);
243 double flengthmax=400;
249 TVector3 distVector = (endPoint-zeroPoint);
252 if(distVector.Mag() < flengthmax)
return false;
253 if(nPoints < 7)
return false;
258 for(
unsigned int j=2;
j<nPoints+1;
j++){
260 distVector = (endPoint-backPoint);
261 if(distVector.Mag()>flengthmax){
270 std::vector<double> msAngles;
272 for(
unsigned int i=skip;
i<(nPoints-(2*fvecsize));
i++){
277 TVector3 vectorA = (bPoint-aPoint);
278 TVector3 vectorB = (cPoint-bPoint);
283 double scalarAB = (vectorA * vectorB);
284 double angleCos = (scalarAB)/(vectorA.Mag()*vectorB.Mag());
285 double angle =
acos (angleCos) * 180.0/TMath::Pi();
287 msAngles.push_back(angle);
288 distVector = (startPoint-cPoint);
289 if(distVector.Mag() >= flengthmax){
290 TVector3 vectorC = (nextPoint-startPoint);
291 TVector3 vectorD = (cPoint-bPoint);
295 double scalarCD = (vectorC * vectorD);
296 double angleExtCos = (scalarCD)/(vectorC.Mag()*vectorD.Mag());
297 angleExt =
acos (angleExtCos) * 180.0/TMath::Pi();
307 for(
unsigned int j=0;
j<msAngles.size();
j++){
308 angleSum += msAngles[
j];
309 angleSumQ += (msAngles[
j]*msAngles[
j]);
310 if(msAngles[
j]>angleMax) angleMax = msAngles[
j];
313 double sizeMS=msAngles.size();
314 double angleMean= (angleSum/sizeMS);
315 angleSigma =
sqrt( (angleSumQ-(2 * angleMean * angleSum)+ (sizeMS * angleMean * angleMean))/sizeMS);
327 std::vector< geo::OfflineChan> xhits, yhits, xhits2, yhits2;
330 TVector3 dirfor = track->
Dir().Unit();
331 TVector3 dirbak = -dirfor;
343 diry = track->
Dir().Y();
345 ratio = (track->
NCell()*1.0) / (1.0*ncell);
346 double cosmicvar = (angle*
angle) - (diry*diry);
347 if(fwddist < 30 && bakdist < 30) cont =
false;
352 if(ncell > 50 && ncell <= 135 && cosmicvar < ((1.1/85)*ncell-1.65)) weakangle =
false;
353 if(ncell > 135 && ncell <= 250 && cosmicvar < ((.6/115)*ncell-.6)) weakangle =
false;
354 if(ncell > 250 && ncell <= 500 && cosmicvar < ((.1/250)*ncell+0.6)) weakangle =
false;
355 if(ncell > 500 && cosmicvar < 0.8) weakangle =
false;
360 if(ncell > 30 && ncell <= 135 && ratio > 0.8 && cosmicvar < ((1/100)*ncell-1.0)) strongangle =
false;
361 if(ncell > 135 && ncell <= 250 && ratio > 0.8 && cosmicvar < ((.45/100)*ncell-.2575)) strongangle =
false;
362 if(ncell > 250 && ncell <= 500 && ratio > 0.83 && cosmicvar < ((.1/200)*ncell+0.6825)) strongangle =
false;
363 if(ncell > 500 && ratio > 0.8 && cosmicvar < 0.9) strongangle =
false;
374 float fMinXFid = -745;
375 float fMinYFid = -745;
378 float fMaxXFid = 745;
379 float fMaxYFid = 745;
382 int stot[6] = {0, 0, 0, 0, 0, 0};
385 for(
unsigned int j = 0;
j < slice->
NCell(); ++
j) {
391 const double x = xyz[0];
392 const double y = xyz[1];
393 const double z = xyz[2];
428 nhitsoutside = *std::max_element(stot,stot+6);
size_t NTrajectoryPoints() const
bool LineIntersection(double x0, double y0, double x1, double y1, double X0, double Y0, double X1, double Y1, double &x, double &y)
Find the intersection between two line segments.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
unsigned short Plane() const
Vertical planes which measure X.
A collection of associated CellHits.
const PlaneGeo * Plane(unsigned int i) const
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
virtual TVector3 Start() const
Horizontal planes which measure Y.
void getBBC(const art::Ptr< rb::Cluster > slice, int &nhitsoutside)
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
unsigned short Cell() const
bool FScatterEstim(art::Ptr< rb::Track > track, bool &fromTrkEnd, double &angleExt, double &angleSigma, double &angleMax, double &angleSum)
double ProjectedLiveDistanceToEdge(TVector3 vertex, TVector3 dir)
void getFits(const art::Ptr< rb::Track > track, double &slope, double &chisq, double &chisqdiff)
double LinFitMinDperp(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Find the best-fit line to a collection of points in 2-D by minimizing the squared perpendicular dista...
double InstrumentedDetFront()
get instrumented detector front of downstream part
olddir
Are we streaming the files via xrootd? txtcmd="cat %s | xargs -n1 samweb2xrootd > xrootd_inFile...
float getScatt(const art::Ptr< rb::Track > track)
double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
In two dimensions, return the perpendicular distance from a point (x0,y0) to a line defined by end po...
float getAngle(TVector3 trackdir)
virtual TVector3 Dir() const
Unit vector describing prong direction.
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Collect Geo headers and supply basic geometry functions.
double InstrumentedDetBack()
get instrumented detector back of downstream part
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
float fitchisq(const std::vector< double > &t, const std::vector< double > &l, int flag, double c)
float fitslope(const std::vector< double > &l, const std::vector< double > &t, const std::vector< double > &q, double &chisq, double &inter, bool linflag)
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
void getCVVars(const art::Ptr< rb::Track > track, int ncell, bool &strongangle, bool &weakangle, bool &cont, float &cvfd, float &cvbd)
Cosmic Rejection PIDs for Numu analysis.
A rawdata::RawDigit with channel information decoded.
float getActivity(art::PtrVector< rb::CellHit > const &trackHits, art::PtrVector< rb::CellHit > const &sliceHits)
assert(nhit_max >=nhit_nbins)
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
TVector3 Stop() const
Position of the final trajectory point.
Encapsulate the geometry of one entire detector (near, far, ndos)