CosRejFxs.cxx
Go to the documentation of this file.
1 #include "CosRej/CosRejFxs.h"
2 
4 #include "Geometry/Geometry.h"
6 #include "GeometryObjects/Geo.h"
8 #include "RecoBase/RecoHit.h"
9 
11 
12 #include "TMath.h"
13 #include "TVector3.h"
14 
15 
16 namespace cosrej
17 {
18 
19  //-------------------------------------------------------
20  CosRejFxs::CosRejFxs() //constructor
21  {
22  }
23 
24  //-------------------------------------------------------
25  CosRejFxs::~CosRejFxs() //destructor
26  {
27  }
28 
29  //-------------------------------------------------------
30  // getAngle: function to get the angle of the lepton wrt the incoming neutrino.
31  // currently the incoming neutrino direction estimate is hardcoded.
32 
33  float CosRejFxs::getAngle(TVector3 trackdir)
34  {
36  TVector3 beamdir = geom->NuMIBeamDirection();
37  float ret = trackdir*beamdir;
38  return ret;
39  }
40 
41  //------------------------------------------------------------------------------
42  // old timing fits - was useful, may re-implement in future
43 
44  void CosRejFxs::getFits(const art::Ptr<rb::Track> track, double &slope, double &chisq, double &chisqdiff)
45  {
46  double xyz[3], getx, gety, getz, intercept;
47  std::vector<double> zs, ts, qs; // vectors for fit: distance along track, time, weight
48  TVector3 tdir = track->Dir();
49  TVector3 tstart = track->Start();
50  TVector3 tend = track->Stop();
51 
54 
55  for(unsigned int hitIdx = 0; hitIdx < track->NCell(); ++hitIdx){
56 
57  art::Ptr<rb::CellHit> chit = track->Cell(hitIdx);
58  geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
59  rb::RecoHit rhit = track->RecoHit(hitIdx);
60 
61  if(!chit->GoodTiming()) {
62  continue;
63  }
64  if(!rhit.IsCalibrated()) continue;
65 
66  if(chit->View()==geo::kX) { // each view is different
67 
68  double px, pz = 1; // make a direction vector perpindicular to track in 2d
69  px = -tdir.Z()/tdir.X();
70  double norm = sqrt(px*px+pz*pz);
71  px = px/norm;
72  pz = pz/norm;
73 
74  geo::LineIntersection(tstart.X(), tstart.Z(), tend.X(), tend.Z(), xyz[0], xyz[2],
75  xyz[0]+50*px, xyz[2]+50*pz, getx, getz); // find intersection with track
76  zs.push_back((getz-tstart.Z())/tdir.Z()); // how far along the track is in 3d
77  }
78 
79  if(chit->View()==geo::kY) {
80 
81  double py, pz = 1; // make a direction vector perindicular to track in 2d
82  py = -tdir.Z()/tdir.Y();
83  double norm = sqrt(py*py+pz*pz);
84  py = py/norm;
85  pz = pz/norm;
86 
87  geo::LineIntersection(tstart.Y(), tstart.Z(), tend.Y(), tend.Z(), xyz[1], xyz[2],
88  xyz[1]+50*py, xyz[2]+50*pz, gety, getz); // find intersection with track
89  zs.push_back((getz-tstart.Z())/tdir.Z());
90  }
91 
92  ts.push_back(rhit.T());
93 
94  if(rhit.IsCalibrated()) {
95  double par1 = -0.377*rhit.PECorr();
96  double par2 = 0.000772*rhit.PECorr()*rhit.PECorr();
97  double par3 = -7.2e-7*rhit.PECorr()*rhit.PECorr()*rhit.PECorr();
98  double par4 = 2.486e-10*rhit.PECorr()*rhit.PECorr()*rhit.PECorr()*rhit.PECorr();
99  double error = 87.3+par1+par2+par3+par4; //timing resolution (ns), fit from timing res plot from Evan
100  if(error<1) error = 1; //keep 0 < 1/error^2 < 1
101  qs.push_back(1/(error*error)); // weight by timing resolution;
102  }
103  else{
104  qs.push_back(1);
105  }
106  }
107 
108  if(zs.size()>1) {
109  slope = fitslope(zs, ts, qs, chisq, intercept, 1);
110  chisqdiff = fitchisq(zs, ts, 0, 30.0) - fitchisq(zs, ts, 1, 30.0);
111  }
112  else {
113  slope = -5;
114  chisqdiff = -5;
115  chisq = -5;
116  }
117 
118  }
119 
120  //---------------------------------------------------------------------------
121 
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) {
123 
124  assert(t.size()==l.size());
125  assert(t.size()>=2);
126 
127  double ztStart[2] = {0.};
128  double ztEnd[2] = {0.};
129  if(linflag==false) {
130  chisq = geo::LinFit(t, l, q,
131  &ztStart[0], &ztStart[1],
132  &ztEnd[0], &ztEnd[1]);
133  }
134  else if(linflag==true) {
135  chisq = geo::LinFitMinDperp(t, l, q,
136  &ztStart[0], &ztStart[1],
137  &ztEnd[0], &ztEnd[1]);
138  }
139  const double ztSlope = (ztEnd[1]-ztStart[1])/(ztEnd[0]-ztStart[0]);
140 
141  inter = ztStart[1]-(ztSlope*ztStart[0]);
142 
143  return ztSlope;
144  }
145 
146  //----------------------------------------------------------------------------------------
147 
148  float CosRejFxs::fitchisq(const std::vector< double > &t, const std::vector< double > &l, int flag, double c) {
149 
150  assert(t.size()==l.size());
151  assert(t.size()>=2);
152 
153  double chi2 = 0.0;
154  double tave = 0;
155  double lave = 0;
156 
157  for (unsigned int i=1; i<t.size(); ++i) {
158  tave += t[i];
159  lave += l[i];
160  }
161 
162  double norm = t.size();
163 
164  tave = tave / norm;
165  lave = lave / norm;
166 
167  double tstart = tave - 500;
168  double tend = tave + 500;
169  double lstart, lend;
170 
171  assert(flag == 0 || flag == 1);
172  if(flag==0) {
173  lstart = lave - (c * 500);
174  lend = lave + (c * 500);
175  }
176  else {
177  lstart = lave + (c * 500);
178  lend = lave - (c * 500);
179  }
180 
181  for (unsigned int i=0; i<t.size(); ++i) {
182  chi2 += geo::DsqrToLine(t[i],l[i],tstart,lstart,tend,lend);
183  }
184  return chi2;
185  }
186 
187  //----------------------------------------------------------------------------------
188 
190  double tscattmax=0, tscatttot=0, tscatt;
191  TVector3 traj1, traj2, dir, olddir;
192 
193  for(unsigned int p=0;p<track->NTrajectoryPoints()-1;p++) {
194  traj1 = track->TrajectoryPoint(p);
195  traj2 = track->TrajectoryPoint(p+1);
196  dir = (traj2-traj1).Unit();
197  if((dir-olddir).Mag()<1e-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;
201  tscatttot += tscatt;
202  }
203  olddir = dir;
204  }
205  return tscatttot; // calculated both maximum and total scatter, total seems more useful, ignore max for now
206  }
207 
208  //-----------------------------------------------------------------------------------
209 
211  // largely copied from NumuEAlg->MakeVertexCluster
212  rb::Cluster vertexCluster;
213  for (unsigned int sliceHit = 0; sliceHit < sliceHits.size(); ++sliceHit){
214  bool match = 0;
215  for (unsigned int trackHit = 0; trackHit < trackHits.size(); ++trackHit){
216  match = (sliceHits[sliceHit] == trackHits[trackHit]);
217  if (match) break;
218  }//End of loop over hits in track
219  if (!match) vertexCluster.Add(sliceHits[sliceHit]);
220  }//End of loop over hits in the slice
221  // std::cout << "New activity clusters: " << trackHits.size() << " " << sliceHits.size() << " " << vertexCluster.NCell();
222  if(vertexCluster.NCell()>0) return vertexCluster.TotalGeV();
223  else return 0;
224  }
225 
226  //--------------------------------------------------------------------------------------
227 
228  bool CosRejFxs::FScatterEstim(art::Ptr<rb::Track> muTrack, bool& fromTrkEnd, double &angleExt, double &angleSigma, double &angleMax, double &angleSum) {
229 
230  /* This fuction make scattering variables for estimating energy.
231 
232  muTrack - muon track
233  fromTrkEnd - specifies which section of the track to use scattering info from
234  0 (default) for the first section of the track
235  1 for last section of the track (NOT an estimate of all energy, only for last piece)
236  output: reco energy for the muon or -1 for tracks too short to use
237 
238  CAVEAT: The reported performance and resolution for this method is ONLY valid for muons with
239  true energy under 4 GeV. This version still returns meaningless values for anything above that!
240  */
241 
242  // Fixed values of flengthmax and fvecsize for this method, could be stored somewhere else
243  double flengthmax=400; // 400 cm of tracks
244  double fvecsize=6; // 6 planes per vector
245 
246  double nPoints = (muTrack->NTrajectoryPoints());
247  TVector3 zeroPoint = (muTrack->TrajectoryPoint(0));
248  TVector3 endPoint = (muTrack->TrajectoryPoint(nPoints-1));
249  TVector3 distVector = (endPoint-zeroPoint);
250 
251  // Check that track is long enough to use
252  if(distVector.Mag() < flengthmax) return false;
253  if(nPoints < 7) return false;
254 
255  int skip = 0;
256  //find point flengthmax cm from end
257  if(fromTrkEnd==1){
258  for(unsigned int j=2; j<nPoints+1; j++){
259  TVector3 backPoint = (muTrack->TrajectoryPoint(nPoints-j));
260  distVector = (endPoint-backPoint);
261  if(distVector.Mag()>flengthmax){
262  skip=nPoints-j;
263  j=nPoints+2;
264  }
265  }
266  }
267 
268  TVector3 startPoint = (muTrack->TrajectoryPoint(skip));
269  TVector3 nextPoint = (muTrack->TrajectoryPoint(skip+fvecsize));
270  std::vector<double> msAngles;
271 
272  for(unsigned int i=skip; i<(nPoints-(2*fvecsize)); i++){
273 
274  TVector3 aPoint = (muTrack->TrajectoryPoint(i));
275  TVector3 bPoint = (muTrack->TrajectoryPoint(i+fvecsize));
276  TVector3 cPoint = (muTrack->TrajectoryPoint(i+(2*fvecsize)));
277  TVector3 vectorA = (bPoint-aPoint);
278  TVector3 vectorB = (cPoint-bPoint);
279 
280  // TVector3 vectorAB = (vectorA * vectorB); // scalar product is not a vector
281  // double angleCos = (vectorAB[0]+vectorAB[1]+vectorAB[2])/(vectorA.Mag()*vectorB.Mag());
282 
283  double scalarAB = (vectorA * vectorB);
284  double angleCos = (scalarAB)/(vectorA.Mag()*vectorB.Mag());
285  double angle = acos (angleCos) * 180.0/TMath::Pi();
286 
287  msAngles.push_back(angle);
288  distVector = (startPoint-cPoint);
289  if(distVector.Mag() >= flengthmax){
290  TVector3 vectorC = (nextPoint-startPoint);
291  TVector3 vectorD = (cPoint-bPoint);
292  // TVector3 vectorCD = (vectorC * vectorD);
293  // double angleExtCos = (vectorCD[0]+vectorCD[1]+vectorCD[2])/(vectorC.Mag()*vectorD.Mag());
294  // angleExt = acos (angleExtCos) * 180.0/TMath::Pi();
295  double scalarCD = (vectorC * vectorD);
296  double angleExtCos = (scalarCD)/(vectorC.Mag()*vectorD.Mag());
297  angleExt = acos (angleExtCos) * 180.0/TMath::Pi();
298 
299  i=nPoints;
300  }
301  }
302 
303  angleSum=0;
304  angleMax=0;
305  double angleSumQ=0;
306  // calculate the rest of the scattering variables
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];
311  }
312 
313  double sizeMS=msAngles.size();
314  double angleMean= (angleSum/sizeMS);
315  angleSigma = sqrt( (angleSumQ-(2 * angleMean * angleSum)+ (sizeMS * angleMean * angleMean))/sizeMS);
316  return true;
317  }
318 
319  //------------------------------------------------------------------------------------------------
320 
321  void CosRejFxs::getCVVars(const art::Ptr<rb::Track> track, int ncell, bool &strongangle, bool &weakangle, bool &cont, float &cvfd, float &cvbd) {
322 
323  cont = true;
324  weakangle = true;
325  strongangle = true;
326 
327  std::vector< geo::OfflineChan> xhits, yhits, xhits2, yhits2;
328  TVector3 start = track->Start();
329  TVector3 end = track->Stop();
330  TVector3 dirfor = track->Dir().Unit();
331  TVector3 dirbak = -dirfor;
332  TVector3 temp;
333  double fwddist, bakdist, diry, ratio, angle;
336 
337  fwddist = livegeom->ProjectedLiveDistanceToEdge(end,dirfor);
338  bakdist = livegeom->ProjectedLiveDistanceToEdge(start,dirbak);
339 
340  cvfd = fwddist;
341  cvbd = bakdist;
342 
343  diry = track->Dir().Y();
344  angle = getAngle(track->Dir());
345  ratio = (track->NCell()*1.0) / (1.0*ncell);
346  double cosmicvar = (angle*angle) - (diry*diry);
347  if(fwddist < 30 && bakdist < 30) cont = false;
348 
349  // weak cut:
350 
351  if(ratio > 0.88) {
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;
356  }
357 
358  // strong cut:
359 
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;
364 
365  } // end GetCVVars fxn
366 
367  //------------------------------------------------------------------------------------------------
368 
369  void CosRejFxs::getBBC(const art::Ptr<rb::Cluster> slice, int &nhitsoutside) {
370 
373 
374  float fMinXFid = -745;
375  float fMinYFid = -745;
376  float fMinZFid = livegeom->InstrumentedDetFront()+24;
377 
378  float fMaxXFid = 745;
379  float fMaxYFid = 745;
380  float fMaxZFid = livegeom->InstrumentedDetBack()-24;
381 
382  int stot[6] = {0, 0, 0, 0, 0, 0};
383 
384  //loop over all hits
385  for(unsigned int j = 0; j < slice->NCell(); ++j) {
386 
387  const rb::CellHit* h = slice->Cell(j).get();
388 
389  double xyz[3];
390  geom->Plane(h->Plane())->Cell(h->Cell())->GetCenter(xyz);
391  const double x = xyz[0];
392  const double y = xyz[1];
393  const double z = xyz[2];
394 
395  // Total up number of cells > Max (X/Y/Z) and < Min (X/Y/Z)
396  //TOP
397  if(y > fMaxYFid) {
398  stot[0]++;
399  continue;
400  }
401  //BOTTOM
402  if(y < fMinYFid) {
403  stot[1]++;
404  continue;
405  }
406  //LEFT
407  if(x > fMaxXFid) {
408  stot[2]++;
409  continue;
410  }
411  //RIGHT
412  if(x < fMinXFid) {
413  stot[3]++;
414  continue;
415  }
416  //FRONT
417  if(z < fMinZFid) {
418  stot[4]++;
419  continue;
420  }
421  //BACK
422  if(z > fMaxZFid) {
423  stot[5]++;
424  continue;
425  }
426  }
427 
428  nhitsoutside = *std::max_element(stot,stot+6);
429  } // end getBBC
430 
431 
432 } //end namespace
float T() const
Definition: RecoHit.h:39
size_t NTrajectoryPoints() const
Definition: Track.h:83
Double_t angle
Definition: plot.C:86
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.
Definition: Geo.cxx:184
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
nPoints
Definition: plotROC.py:88
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...
Definition: Geo.cxx:373
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1 * ratio(TH1 *h1, TH1 *h2)
Vertical planes which measure X.
Definition: PlaneGeo.h:28
float Z() const
Definition: RecoHit.h:38
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
T acos(T number)
Definition: d0nt_math.hpp:54
double chi2()
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void getBBC(const art::Ptr< rb::Cluster > slice, int &nhitsoutside)
Definition: CosRejFxs.cxx:369
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
unsigned short Cell() const
Definition: CellHit.h:40
bool FScatterEstim(art::Ptr< rb::Track > track, bool &fromTrkEnd, double &angleExt, double &angleSigma, double &angleMax, double &angleSum)
Definition: CosRejFxs.cxx:228
double ProjectedLiveDistanceToEdge(TVector3 vertex, TVector3 dir)
void getFits(const art::Ptr< rb::Track > track, double &slope, double &chisq, double &chisqdiff)
Definition: CosRejFxs.cxx:44
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...
Definition: Geo.cxx:449
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...
Definition: runNovaSAM.py:800
float getScatt(const art::Ptr< rb::Track > track)
Definition: CosRejFxs.cxx:189
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...
Definition: Geo.cxx:338
float getAngle(TVector3 trackdir)
Definition: CosRejFxs.cxx:33
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
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.
Definition: RecoHit.cxx:35
const double j
Definition: BetheBloch.cxx:29
TVector3 Unit() const
z
Definition: test.py:28
float fitchisq(const std::vector< double > &t, const std::vector< double > &l, int flag, double c)
Definition: CosRejFxs.cxx:148
size_type size() const
Definition: PtrVector.h:308
float fitslope(const std::vector< double > &l, const std::vector< double > &t, const std::vector< double > &q, double &chisq, double &inter, bool linflag)
Definition: CosRejFxs.cxx:122
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.
Definition: Cluster.cxx:145
void getCVVars(const art::Ptr< rb::Track > track, int ncell, bool &strongangle, bool &weakangle, bool &cont, float &cvfd, float &cvbd)
Definition: CosRejFxs.cxx:321
Cosmic Rejection PIDs for Numu analysis.
Definition: FillParentInfo.h:9
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
Float_t norm
T const * get() const
Definition: Ptr.h:321
TDirectory * dir
Definition: macro.C:5
void geom(int which=0)
Definition: geom.C:163
float getActivity(art::PtrVector< rb::CellHit > const &trackHits, art::PtrVector< rb::CellHit > const &sliceHits)
Definition: CosRejFxs.cxx:210
assert(nhit_max >=nhit_nbins)
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
float PECorr() const
Definition: RecoHit.cxx:47
bool GoodTiming() const
Definition: CellHit.h:48
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Float_t e
Definition: plot.C:35
virtual ~CosRejFxs()
Definition: CosRejFxs.cxx:25
Encapsulate the geometry of one entire detector (near, far, ndos)