CosRejFxs.cxx

Go to the documentation of this file.

#include "CosRej/CosRejFxs.h"

#include "Calibrator/Calibrator.h"
#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/RecoHit.h"

#include "art/Framework/Services/Registry/ServiceHandle.h"

#include "TMath.h"
#include "TVector3.h"


namespace cosrej
{

  //-------------------------------------------------------
  CosRejFxs::CosRejFxs() //constructor
  {
  }

  //-------------------------------------------------------
  CosRejFxs::~CosRejFxs() //destructor
  {
  }

  //-------------------------------------------------------
  // getAngle: function to get the angle of the lepton wrt the incoming neutrino.
  // currently the incoming neutrino direction estimate is hardcoded.

  float CosRejFxs::getAngle(TVector3 trackdir)
  {
    art::ServiceHandle<geo::Geometry> geom;
    TVector3 beamdir = geom->NuMIBeamDirection();
    float ret = trackdir*beamdir;
    return ret;
  }

  //------------------------------------------------------------------------------
  // old timing fits - was useful, may re-implement in future

  void CosRejFxs::getFits(const art::Ptr<rb::Track> track, double &slope, double &chisq, double &chisqdiff)
  {
    double xyz[3], getx, gety, getz, intercept;
    std::vector<double> zs, ts, qs; // vectors for fit: distance along track, time, weight
    TVector3 tdir = track->Dir();
    TVector3 tstart = track->Start();
    TVector3 tend = track->Stop();

    art::ServiceHandle<geo::Geometry>  geom;
    art::ServiceHandle<calib::Calibrator> cal;

    for(unsigned int hitIdx = 0; hitIdx < track->NCell(); ++hitIdx){

      art::Ptr<rb::CellHit> chit = track->Cell(hitIdx);
      geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
      rb::RecoHit rhit = track->RecoHit(hitIdx);

      if(!chit->GoodTiming()) {
        continue;
      }
      if(!rhit.IsCalibrated()) continue;

      if(chit->View()==geo::kX) { // each view is different

        double px, pz = 1;  // make a direction vector perpindicular to track in 2d
        px = -tdir.Z()/tdir.X();
        double norm = sqrt(px*px+pz*pz);
        px = px/norm;
        pz = pz/norm;

        geo::LineIntersection(tstart.X(), tstart.Z(), tend.X(), tend.Z(), xyz[0], xyz[2],
                              xyz[0]+50*px, xyz[2]+50*pz, getx, getz);  // find intersection with track
        zs.push_back((getz-tstart.Z())/tdir.Z());  // how far along the track is in 3d
      }

      if(chit->View()==geo::kY) {

        double py, pz = 1;  // make a direction vector perindicular to track in 2d
        py = -tdir.Z()/tdir.Y();
        double norm = sqrt(py*py+pz*pz);
        py = py/norm;
        pz = pz/norm;

        geo::LineIntersection(tstart.Y(), tstart.Z(), tend.Y(), tend.Z(), xyz[1], xyz[2],
                              xyz[1]+50*py, xyz[2]+50*pz, gety, getz);  // find intersection with track
        zs.push_back((getz-tstart.Z())/tdir.Z());
      }

      ts.push_back(rhit.T());

      if(rhit.IsCalibrated()) {
        double par1 = -0.377*rhit.PECorr();
        double par2 = 0.000772*rhit.PECorr()*rhit.PECorr();
        double par3 = -7.2e-7*rhit.PECorr()*rhit.PECorr()*rhit.PECorr();
        double par4 = 2.486e-10*rhit.PECorr()*rhit.PECorr()*rhit.PECorr()*rhit.PECorr();
        double error = 87.3+par1+par2+par3+par4;  //timing resolution (ns), fit from timing res plot from Evan
        if(error<1) error = 1; //keep 0 < 1/error^2 < 1
        qs.push_back(1/(error*error));  // weight by timing resolution;
      }
      else{
        qs.push_back(1);
      }
    }

    if(zs.size()>1) {
      slope = fitslope(zs, ts, qs, chisq, intercept, 1);
      chisqdiff = fitchisq(zs, ts, 0, 30.0) - fitchisq(zs, ts, 1, 30.0);
    }
    else {
      slope = -5;
      chisqdiff = -5;
      chisq = -5;
    }

  }

  //---------------------------------------------------------------------------

  float CosRejFxs::fitslope(const std::vector< double > &l, const std::vector< double > &t, const std::vector< double > &q, double &chisq, double &inter, bool linflag) {

    assert(t.size()==l.size());
    assert(t.size()>=2);

    double ztStart[2] = {0.};
    double ztEnd[2]   = {0.};
    if(linflag==false) {
      chisq = geo::LinFit(t, l, q,
                          &ztStart[0], &ztStart[1],
                          &ztEnd[0], &ztEnd[1]);
    }
    else if(linflag==true) {
      chisq = geo::LinFitMinDperp(t, l, q,
                                  &ztStart[0], &ztStart[1],
                                  &ztEnd[0], &ztEnd[1]);
    }
    const double ztSlope = (ztEnd[1]-ztStart[1])/(ztEnd[0]-ztStart[0]);

    inter = ztStart[1]-(ztSlope*ztStart[0]);

    return ztSlope;
  }

  //----------------------------------------------------------------------------------------

  float CosRejFxs::fitchisq(const std::vector< double > &t, const std::vector< double > &l, int flag, double c) {

    assert(t.size()==l.size());
    assert(t.size()>=2);

    double chi2 = 0.0;
    double tave = 0;
    double lave = 0;

    for (unsigned int i=1; i<t.size(); ++i) {
      tave += t[i];
      lave += l[i];
    }

    double norm = t.size();

    tave = tave / norm;
    lave = lave / norm;

    double tstart = tave - 500;
    double tend = tave + 500;
    double lstart, lend;

    assert(flag == 0 || flag == 1);
    if(flag==0) {
      lstart = lave - (c * 500);
      lend = lave + (c * 500);
    }
    else {
      lstart = lave + (c * 500);
      lend = lave - (c * 500);
    }

    for (unsigned int i=0; i<t.size(); ++i) {
      chi2 += geo::DsqrToLine(t[i],l[i],tstart,lstart,tend,lend);
    }
    return chi2;
  }

  //----------------------------------------------------------------------------------

  float CosRejFxs::getScatt(const art::Ptr<rb::Track> track) {
    double tscattmax=0, tscatttot=0, tscatt;
    TVector3 traj1, traj2, dir, olddir;

    for(unsigned int p=0;p<track->NTrajectoryPoints()-1;p++) {
      traj1 = track->TrajectoryPoint(p);
      traj2 = track->TrajectoryPoint(p+1);
      dir = (traj2-traj1).Unit();
      if((dir-olddir).Mag()<1e-7) continue;
      if(p>0&&dir.Mag()<=1.0&&olddir.Mag()<=1.0) {
        tscatt = std::acos(dir.Dot(olddir)) * 180 / TMath::Pi();
        if(tscatt > tscattmax) tscattmax = tscatt;
        tscatttot += tscatt;
      }
      olddir = dir;
    }
    return tscatttot;  // calculated both maximum and total scatter, total seems more useful, ignore max for now
  }

  //-----------------------------------------------------------------------------------

  float CosRejFxs::getActivity(art::PtrVector< rb::CellHit > const& trackHits, art::PtrVector< rb::CellHit > const& sliceHits) {
    // largely copied from NumuEAlg->MakeVertexCluster
    rb::Cluster vertexCluster;
    for (unsigned int sliceHit = 0; sliceHit < sliceHits.size(); ++sliceHit){
      bool match = 0;
      for (unsigned int trackHit = 0; trackHit < trackHits.size(); ++trackHit){
        match = (sliceHits[sliceHit] == trackHits[trackHit]);
        if (match) break;
      }//End of loop over hits in track
      if (!match) vertexCluster.Add(sliceHits[sliceHit]);
    }//End of loop over hits in the slice
    //    std::cout << "New activity clusters: " << trackHits.size() << " " << sliceHits.size() << " " << vertexCluster.NCell();
    if(vertexCluster.NCell()>0) return vertexCluster.TotalGeV();
    else return 0;
  }

  //--------------------------------------------------------------------------------------

  bool CosRejFxs::FScatterEstim(art::Ptr<rb::Track> muTrack, bool& fromTrkEnd, double &angleExt, double &angleSigma, double &angleMax, double &angleSum) {

    /*  This fuction make scattering variables for estimating energy.

      muTrack - muon track
      fromTrkEnd - specifies which section of the track to use scattering info from
        0 (default) for the first section of the track
        1 for last section of the track (NOT an estimate of all energy, only for last piece)
      output: reco energy for the muon or -1 for tracks too short to use

      CAVEAT: The reported performance and resolution for this method is ONLY valid for muons with
      true energy under 4 GeV. This version still returns meaningless values for anything above that!
    */

    // Fixed values of flengthmax and fvecsize for this method, could be stored somewhere else
    double flengthmax=400;          // 400 cm of tracks
    double fvecsize=6;              // 6 planes per vector

    double nPoints      = (muTrack->NTrajectoryPoints());
    TVector3 zeroPoint  = (muTrack->TrajectoryPoint(0));
    TVector3 endPoint   = (muTrack->TrajectoryPoint(nPoints-1));
    TVector3 distVector = (endPoint-zeroPoint);

    // Check that track is long enough to use
    if(distVector.Mag() < flengthmax) return false;
    if(nPoints < 7) return false;

    int skip = 0;
    //find point flengthmax cm from end
    if(fromTrkEnd==1){
      for(unsigned int j=2; j<nPoints+1; j++){
        TVector3 backPoint = (muTrack->TrajectoryPoint(nPoints-j));
        distVector = (endPoint-backPoint);
        if(distVector.Mag()>flengthmax){
          skip=nPoints-j;
          j=nPoints+2;
        }
      }
    }

    TVector3 startPoint = (muTrack->TrajectoryPoint(skip));
    TVector3 nextPoint  = (muTrack->TrajectoryPoint(skip+fvecsize));
    std::vector<double> msAngles;

    for(unsigned int i=skip; i<(nPoints-(2*fvecsize)); i++){

      TVector3 aPoint   = (muTrack->TrajectoryPoint(i));
      TVector3 bPoint   = (muTrack->TrajectoryPoint(i+fvecsize));
      TVector3 cPoint   = (muTrack->TrajectoryPoint(i+(2*fvecsize)));
      TVector3 vectorA  = (bPoint-aPoint);
      TVector3 vectorB  = (cPoint-bPoint);

      //      TVector3 vectorAB = (vectorA * vectorB); // scalar product is not a vector
      //      double angleCos   = (vectorAB[0]+vectorAB[1]+vectorAB[2])/(vectorA.Mag()*vectorB.Mag());

      double scalarAB = (vectorA * vectorB);
      double angleCos   = (scalarAB)/(vectorA.Mag()*vectorB.Mag());
      double angle      = acos (angleCos) * 180.0/TMath::Pi();

      msAngles.push_back(angle);
      distVector = (startPoint-cPoint);
      if(distVector.Mag() >= flengthmax){
        TVector3 vectorC   = (nextPoint-startPoint);
        TVector3 vectorD   = (cPoint-bPoint);
        // TVector3 vectorCD  = (vectorC * vectorD);
        // double angleExtCos = (vectorCD[0]+vectorCD[1]+vectorCD[2])/(vectorC.Mag()*vectorD.Mag());
        // angleExt = acos (angleExtCos) * 180.0/TMath::Pi();
        double scalarCD  = (vectorC * vectorD);
        double angleExtCos = (scalarCD)/(vectorC.Mag()*vectorD.Mag());
        angleExt = acos (angleExtCos) * 180.0/TMath::Pi();

        i=nPoints;
      }
    }

    angleSum=0;
    angleMax=0;
    double angleSumQ=0;
    // calculate the rest of the scattering variables
    for(unsigned int j=0; j<msAngles.size(); j++){
      angleSum += msAngles[j];
      angleSumQ += (msAngles[j]*msAngles[j]);
      if(msAngles[j]>angleMax) angleMax = msAngles[j];
    }

    double sizeMS=msAngles.size();
    double angleMean= (angleSum/sizeMS);
    angleSigma = sqrt( (angleSumQ-(2 * angleMean * angleSum)+ (sizeMS * angleMean * angleMean))/sizeMS);
    return true;
  }

  //------------------------------------------------------------------------------------------------

  void CosRejFxs::getCVVars(const art::Ptr<rb::Track> track, int ncell, bool &strongangle, bool &weakangle, bool &cont, float &cvfd, float &cvbd) {

    cont = true;
    weakangle = true;
    strongangle = true;

    std::vector< geo::OfflineChan> xhits, yhits, xhits2, yhits2;
    TVector3 start = track->Start();
    TVector3 end = track->Stop();
    TVector3 dirfor = track->Dir().Unit();
    TVector3 dirbak = -dirfor;
    TVector3 temp;
    double fwddist, bakdist, diry, ratio, angle;
    art::ServiceHandle<geo::Geometry>  geom;
    art::ServiceHandle<geo::LiveGeometry>  livegeom;

    fwddist = livegeom->ProjectedLiveDistanceToEdge(end,dirfor);
    bakdist = livegeom->ProjectedLiveDistanceToEdge(start,dirbak);

    cvfd = fwddist;
    cvbd = bakdist;

    diry = track->Dir().Y();
    angle = getAngle(track->Dir());
    ratio = (track->NCell()*1.0) / (1.0*ncell);
    double cosmicvar = (angle*angle) - (diry*diry);
    if(fwddist < 30 && bakdist < 30) cont = false;

    // weak cut:

    if(ratio > 0.88) {
      if(ncell > 50  && ncell <= 135 && cosmicvar < ((1.1/85)*ncell-1.65)) weakangle = false;
      if(ncell > 135 && ncell <= 250 && cosmicvar < ((.6/115)*ncell-.6)) weakangle = false;
      if(ncell > 250 && ncell <= 500 && cosmicvar < ((.1/250)*ncell+0.6)) weakangle = false;
      if(ncell > 500 && cosmicvar < 0.8) weakangle = false;
    }

    // strong cut:

    if(ncell > 30 && ncell <= 135 && ratio > 0.8 && cosmicvar < ((1/100)*ncell-1.0)) strongangle = false;
    if(ncell > 135 && ncell <= 250 && ratio > 0.8 && cosmicvar < ((.45/100)*ncell-.2575)) strongangle = false;
    if(ncell > 250 && ncell <= 500 && ratio > 0.83 && cosmicvar < ((.1/200)*ncell+0.6825)) strongangle = false;
    if(ncell > 500 && ratio > 0.8 && cosmicvar < 0.9) strongangle = false;

  } // end GetCVVars fxn

  //------------------------------------------------------------------------------------------------

  void CosRejFxs::getBBC(const art::Ptr<rb::Cluster> slice, int &nhitsoutside) {

    art::ServiceHandle<geo::Geometry>  geom;
    art::ServiceHandle<geo::LiveGeometry>  livegeom;

    float fMinXFid = -745;
    float fMinYFid = -745;
    float fMinZFid = livegeom->InstrumentedDetFront()+24;

    float fMaxXFid =  745;
    float fMaxYFid =  745;
    float fMaxZFid = livegeom->InstrumentedDetBack()-24;

    int stot[6] = {0, 0, 0, 0, 0, 0};

    //loop over all hits
    for(unsigned int j = 0; j < slice->NCell(); ++j) {

      const rb::CellHit* h = slice->Cell(j).get();

      double xyz[3];
      geom->Plane(h->Plane())->Cell(h->Cell())->GetCenter(xyz);
      const double x = xyz[0];
      const double y = xyz[1];
      const double z = xyz[2];

      // Total up number of cells > Max (X/Y/Z) and < Min (X/Y/Z)
      //TOP
      if(y > fMaxYFid) {
        stot[0]++;
        continue;
      }
      //BOTTOM
      if(y < fMinYFid) {
        stot[1]++;
        continue;
      }
      //LEFT
      if(x > fMaxXFid) {
        stot[2]++;
        continue;
      }
      //RIGHT
      if(x < fMinXFid) {
        stot[3]++;
        continue;
      }
      //FRONT
      if(z < fMinZFid) {
        stot[4]++;
        continue;
      }
      //BACK
      if(z > fMaxZFid) {
        stot[5]++;
        continue;
      }
    }

    nhitsoutside = *std::max_element(stot,stot+6);
  } // end getBBC


} //end namespace