FMMTracker_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       FMMTracker
// Module Type: producer
// File:        FMMTracker_module.cc
//
// Generated at Wed Dec 18 15:41:32 2013 by Zukai Wang using artmod
// from cetpkgsupport v1_04_02.
////////////////////////////////////////////////////////////////////////
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"

#include "MCCheater/BackTracker.h"
#include "Calibrator/Calibrator.h"
#include "Simulation/Particle.h"
#include "Simulation/ParticleNavigator.h"
#include "art/Framework/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "Utilities/func/MathUtil.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Services/Optional/TFileService.h"

#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/PlaneGeo.h"

#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"

#include "TVector3.h"
#include "TTree.h"

#include <algorithm> 
#include <utility>
#include <cmath>
#include <iostream>

namespace zcl {
  class FMMTracker;
}

class zcl::FMMTracker : public art::EDProducer {
public:
  explicit FMMTracker(fhicl::ParameterSet const & pset);
  virtual ~FMMTracker();

  void produce(art::Event & evt);

  bool Tracking2D(art::ServiceHandle<geo::Geometry> const & geom,
                  const art::PtrVector<rb::CellHit> & hitlist, 
                  int view,
                  rb::Track & track,
                  std::vector<TVector3> & poslist);  

  double recoT(TVector3 position, 
               std::vector<TVector3>& poslist, 
               std::vector<double>& PEs,
               std::vector<double>& tlist);

  double calcW(double z, double w1, double w2, double z1, double z2); 

  private:
  std::string fClusterInput;    ///< Input folder from cluster reco
  double _dSqrMax;      //Determine whether this hit will be included in the track
  int _PEcut;         //whether this hit is qualified to place a weight in the timing of start or stop point
  unsigned _minHits;
  //float _deltaT;
  //float _deltaZ;

};

zcl::FMMTracker::FMMTracker(fhicl::ParameterSet const & p) :
  fClusterInput      (p.get< std::string > ("ClusterInput")      ),
  _dSqrMax           (200),
  _PEcut             (500),
  _minHits           (3)//,     //At least 3 hits on each view
  // _deltaT            (1500),  //ns
  // _deltaZ            (30)     //cm
{
  // Call appropriate Produces<>() functions here.
  produces< std::vector<rb::Track> >();
}

zcl::FMMTracker::~FMMTracker()
{
  // Clean up dynamic memory and other resources here.
}

void zcl::FMMTracker::produce(art::Event & evt)
{
  art::ServiceHandle<geo::Geometry> geom;
  art::ServiceHandle<calib::Calibrator> calib;

  //Now, start reco tracks!
  art::Handle< std::vector<rb::Cluster> > slices;
  evt.getByLabel(fClusterInput, slices);
  std::unique_ptr< std::vector<rb::Track> > trackcol(new std::vector<rb::Track>);

  //Instead of splitting the hits into 2 views, just directly take the entire 3D slice:
  for (unsigned sliceIdx = 0; sliceIdx != slices->size(); ++sliceIdx) {
    const rb::Cluster& slice = (*slices)[sliceIdx];

    const art::PtrVector<rb::CellHit>& Xhitlist = slice.XCells();
    //During the stage of the MMSlicer, the hits should already be sorted by plane number
    const art::PtrVector<rb::CellHit>& Yhitlist = slice.YCells();

    int px_min = 999;
    int px_max = -999;
    double tx_min = 9e9;
    double tx_max = -9e9;

    for (unsigned i=0; i!= Xhitlist.size(); ++i) {
      const art::Ptr<rb::CellHit> chit = Xhitlist[i];
      if (tx_min > chit->TNS()) tx_min = chit->TNS();
      if (tx_max < chit->TNS()) tx_max = chit->TNS();
      if (px_min > chit->Plane()) px_min = chit->Plane();
      if (px_max < chit->Plane()) px_max = chit->Plane();
    }

    int py_min = 999;
    int py_max = -999;
    double ty_min = 9e9;
    double ty_max = -9e9;

    for (unsigned i=0; i!= Yhitlist.size(); ++i) {
      const art::Ptr<rb::CellHit> chit = Yhitlist[i];
      if (ty_min > chit->TNS()) ty_min = chit->TNS();
      if (ty_max < chit->TNS()) ty_max = chit->TNS();
      if (py_min > chit->Plane()) py_min = chit->Plane();
      if (py_max < chit->Plane()) py_max = chit->Plane();
    }

    //Matching test:
    if (std::min(tx_max, ty_max) - std::max(tx_min,ty_min) < 0) continue;
    if (std::min(px_max, py_max) - std::max(px_min,py_min) < 3) continue;

    //    rb::SortByPlane(Yhitlist);
    rb::Track xtrack, ytrack;
    //    TVector3 pX1, pX2, pY1, pY2;
    std::vector<TVector3> Xp_list,Yp_list; //X(Y),Z
    std::vector<TVector3> trjList;         //X,Y,Z
    std::vector<double> t_list;

    //Either X or Y track is not reconstructed, just skip to the next slice
    // geom: hitlist: view: track: positions list:  
    if (!Tracking2D(geom,Xhitlist,0,xtrack,Xp_list) ) 
      continue;

    if (!Tracking2D(geom,Yhitlist,1,ytrack,Yp_list) ) 
      continue;

    //Merging: and setting W position for each reco hit, to get the calibrated time and energy
    rb::Track FMMtrack3D;
    //Reconstruct the 3D info of the start and end points

    TVector3 lowZ_P, highZ_P;
    //    double   lowZ_T, highZ_T;

    unsigned nxhits_RC = Xp_list.size();
    unsigned nyhits_RC = Yp_list.size();

    for (unsigned i=0; i!= nxhits_RC; ++i) {
      const art::Ptr<rb::CellHit> chit = xtrack.Cell(i);
 
      rb::RecoHit rhit(calib->MakeRecoHit(*chit,
                                          calcW(Xp_list[i].Y(),
                                                Yp_list[0].X(),
                                                Yp_list[nyhits_RC-1].X(),
                                                Yp_list[0].Y(),
                                                Yp_list[nyhits_RC-1].Y() ) ) );
      if (!rhit.IsCalibrated() )
        continue;

      trjList.emplace_back(TVector3(rhit.X(),rhit.Y(),rhit.Z()));
      t_list.emplace_back(rhit.T());
      FMMtrack3D.AppendTrajectoryPoint(TVector3(rhit.X(),rhit.Y(),rhit.Z()) );

      //     std::cout<<"Xview:  "<<"X: "<<rhit.X()<<" Y:"<<rhit.Y()<<" Z: "<<rhit.Z()<<std::endl;

      FMMtrack3D.Add(chit);   
    }

    //X view finished, predetermine the lowZ and highZ:
    /*
    lowZ_P  = trjList[0];
    lowZ_T  = t_list[0];
    highZ_P = trjList[nxhits_RC-1];
    highZ_T = t_list[nxhits_RC-1];
    */
    //    std::cout<<"Y hits on track: "<<Yp_list.size()<<std::endl;
    for (unsigned i=0; i!= nyhits_RC; ++i) {
      const art::Ptr<rb::CellHit> chit = ytrack.Cell(i);
      rb::RecoHit rhit(calib->MakeRecoHit(*chit,
                                          calcW(Yp_list[i].Y(),
                                                Xp_list[0].X(),
                                                Xp_list[nxhits_RC-1].X(),
                                                Xp_list[0].Y(),
                                                Xp_list[nxhits_RC-1].Y() ) ) );
  
      if (!rhit.IsCalibrated() )
        continue;
      //     std::cout<<"Yview:  "<<"X: "<<rhit.X()<<" Y:"<<rhit.Y()<<" Z: "<<rhit.Z()<<std::endl;

      trjList.emplace_back(TVector3(rhit.X(),rhit.Y(),rhit.Z()));
      t_list.emplace_back(rhit.T());
      FMMtrack3D.AppendTrajectoryPoint(TVector3(rhit.X(),rhit.Y(),rhit.Z()) );
      FMMtrack3D.Add(chit);   
    }
    /*
    if (lowZ_P.Z() > trjList[nxhits_RC].Z()) {
      lowZ_P = trjList[nxhits_RC];
      lowZ_T = t_list[nxhits_RC];
    }
   
    if (highZ_P.Z() < trjList[nxhits_RC+nyhits_RC-1].Z()) {
      highZ_P = trjList[nxhits_RC+nyhits_RC-1];
      highZ_T = t_list[nxhits_RC+nyhits_RC-1];
    }
 
    //Now, determine the start and end time (simultaneously you know which is the start point):
    lowZ_T   = recoT(lowZ_P,trjList,pe_list,t_list);
    highZ_T  = recoT(highZ_P,trjList,pe_list,t_list);
   

    if (lowZ_T < highZ_T) {
      FMMtrack3D.SetStart(lowZ_P);
      FMMtrack3D.SetDir((highZ_P-lowZ_P).Unit());
    }
    else {
      FMMtrack3D.SetStart(highZ_P);
      FMMtrack3D.SetDir((lowZ_P-highZ_P).Unit());
    }
    */
    trackcol->emplace_back(FMMtrack3D); 
    //    std::cout<<"reconstructed 3D tracks: "<<trackcol->size()<<std::endl;
  }
  evt.put(std::move(trackcol));
} //end production function


//Predetermine the first and last hit in this 2D track
bool zcl::FMMTracker::Tracking2D(art::ServiceHandle<geo::Geometry> const & geom,
                                 const art::PtrVector<rb::CellHit> & hitlist,
                                 int view,                        //For index convenience
                                 rb::Track & track,
                                 std::vector<TVector3> & plist    //X(Y), Z, T
                                 ) 
{
  double xyz[3];
  std::vector<double> x,y,w,t;

  if (hitlist.size()<_minHits) return false;

  //First time loop: store xyw
  for (unsigned i=0; i!= hitlist.size(); ++i) {
    const art::Ptr<rb::CellHit> chit = hitlist[i];


    w.emplace_back(chit->ADC() );
    geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
    x.emplace_back(xyz[2]);
    y.emplace_back(xyz[view]);
    t.emplace_back(chit->TNS());
  } 
  //Doing Linfit to determine the start and end position:
  //But actually I just use this to drop non-in-line hits in the cluster 
  double x1,y1,x2,y2,dSq;
  dSq=geo::LinFitMinDperp(x,y,w,&x1,&x2,&y1,&y2);
  
  //Second time loop: drop hits not on the line
  for (unsigned i=0; i!= hitlist.size(); ++i) {
    dSq = geo::DsqrToLine(x[i],y[i],x1,x2,y1,y2);
    if (dSq<_dSqrMax) {
      track.Add(hitlist[i]);
      plist.emplace_back( TVector3(y[i], x[i], t[i]) );
    }
  }

  if (track.NCell()<_minHits) return false;
  else return true;

}

double zcl::FMMTracker::recoT(TVector3 position, 
                              std::vector<TVector3>& poslist, 
                              std::vector<double>& PEs,
                              std::vector<double>& tlist) 
{
  double T  = 0;
  double wt = 0;
  for (unsigned i=0; i!= PEs.size(); ++i) {
    if (PEs[i] < _PEcut) continue;
    double weight = exp(-(poslist[i]-position).Mag2()); 
    T  = weight*tlist[i]; 
    wt += weight;
  }
  T = T/wt;
  return T;
}

double zcl::FMMTracker::calcW(double z, double w1, double w2, double z1, double z2) {
  double w =  (w2*(z-z1)+w1*(z2-z))/(z2-z1);
  if (w > 761) w = 761;
  else if (w < -761) w = -761;

  return w;
}



DEFINE_ART_MODULE(zcl::FMMTracker)