trkf::Track3DKalman Class Reference

Inheritance diagram for trkf::Track3DKalman:

Public Member Functions
	Track3DKalman (fhicl::ParameterSet const &pset)

Private Member Functions
void	produce (art::Event &evt) override
void	beginJob () override
void	endJob () override

Private Attributes
std::string	fSpacePtsModuleLabel
std::string	fGenieGenModuleLabel
std::string	fG4ModuleLabel
bool	fGenfPRINT
TTree *	tree
TMatrixT< Double_t > *	stMCT
TMatrixT< Double_t > *	covMCT
TMatrixT< Double_t > *	stREC
TMatrixT< Double_t > *	covREC
Float_t	chi2
Float_t	chi2ndf
Float_t *	fpRECt3D
Float_t *	fpRECL
Float_t *	fpREC
Float_t *	fpMCT
int	nfail
int	ndf
unsigned int	evtt
unsigned int	nTrks
unsigned int	fptsNo
Float_t *	fshx
Float_t *	fshy
Float_t *	fshz
unsigned int	fDimSize
std::vector< double >	fPosErr
std::vector< double >	fMomErr
std::vector< double >	fMomStart
genf::GFAbsTrackRep *	repMC
genf::GFAbsTrackRep *	rep
CLHEP::HepRandomEngine &	fEngine

Detailed Description

Definition at line 62 of file Track3DKalman_module.cc.

Constructor & Destructor Documentation

trkf::Track3DKalman::Track3DKalman ( fhicl::ParameterSet const &  pset )

explicit

Definition at line 121 of file Track3DKalman_module.cc.

  : EDProducer{pset}
  , fSpacePtsModuleLabel{pset.get< std::string >("SpacePtsModuleLabel")}
  , fGenieGenModuleLabel{pset.get< std::string >("GenieGenModuleLabel")}
  , fG4ModuleLabel{pset.get< std::string >("G4ModuleLabel")}
  , fGenfPRINT{pset.get< bool >("GenfPRINT")}
  , fPosErr{pset.get< std::vector < double >  >("PosErr3")}   // resolution. cm
  , fMomErr{pset.get< std::vector < double >  >("MomErr3")}   // GeV
  , fMomStart{pset.get< std::vector < double >  >("MomStart3")} // Will be unit norm'd.
  // create a default random engine; obtain the random seed from NuRandomService,
  // unless overridden in configuration with key "Seed"
  , fEngine(art::ServiceHandle<rndm::NuRandomService>()->createEngine(*this, pset, "Seed"))
{
  produces< std::vector<recob::Track> >();
  produces< std::vector<recob::SpacePoint>              >();
  produces< art::Assns<recob::Track, recob::Cluster>    >();
  produces< art::Assns<recob::Track, recob::SpacePoint> >();
  produces< art::Assns<recob::Track, recob::Hit>        >();
}

Member Function Documentation

void trkf::Track3DKalman::beginJob ( )

overrideprivate

Definition at line 142 of file Track3DKalman_module.cc.

{
  art::ServiceHandle<art::TFileService const> tfs;

  stMCT  = new TMatrixT<Double_t>(5,1);
  covMCT = new TMatrixT<Double_t>(5,5);
  stREC  = new TMatrixT<Double_t>(5,1);
  covREC = new TMatrixT<Double_t>(5,5);

  fpMCT = new Float_t[4];
  fpREC = new Float_t[4];
  fpRECL = new Float_t[4];
  fpRECt3D = new Float_t[4];
  fDimSize = 400; // if necessary will get this from pset in constructor.

  fshx = new Float_t[fDimSize];
  fshy = new Float_t[fDimSize];
  fshz = new Float_t[fDimSize];

  tree = tfs->make<TTree>("GENFITttree","GENFITttree");
  //tree->Branch("stMCT",&stMCT,"stMCT[5]/F"); // "TMatrixT<Double_t>"

  tree->Branch("stMCT","TMatrixD",&stMCT,64000,0);
  //tree->Branch("covMCT",&covMCT,"covMCT[25]/F");
  tree->Branch("covMCT","TMatrixD",&covMCT,64000,0);
  //tree->Branch("stREC",&stREC,"stREC[5]/F");
  tree->Branch("stREC","TMatrixD",&stREC,64000,0);
  //tree->Branch("covREC",&covREC,"covREC[25]/F");
  tree->Branch("covREC","TMatrixD",&covREC,64000,0);


  tree->Branch("chi2",&chi2,"chi2/F");
  tree->Branch("nfail",&nfail,"nfail/I");
  tree->Branch("ndf",&ndf,"ndf/I");
  tree->Branch("evtNo",&evtt,"evtNo/I");
  tree->Branch("chi2ndf",&chi2ndf,"chi2ndf/F");

  tree->Branch("trkNo",&nTrks,"trkNo/I");
  tree->Branch("ptsNo",&fptsNo,"ptsNo/I");
  tree->Branch("shx",fshx,"shx[ptsNo]/F");
  tree->Branch("shy",fshy,"shy[ptsNo]/F");
  tree->Branch("shz",fshz,"shz[ptsNo]/F");

  tree->Branch("pMCT",fpMCT,"pMCT[4]/F");
  tree->Branch("pRECKalF",fpREC,"pRECKalF[4]/F");
  tree->Branch("pRECKalL",fpRECL,"pRECKalL[4]/F");
  tree->Branch("pRECt3D",fpRECt3D,"pRECt3D[4]/F");
}

void trkf::Track3DKalman::endJob ( )

overrideprivate

Definition at line 192 of file Track3DKalman_module.cc.

{
  if (!rep) delete rep;
  if (!repMC) delete repMC;
}

void trkf::Track3DKalman::produce ( art::Event &  evt )

overrideprivate

Todo:

: remove this statement, there is no place for checks on isRealData in reconstruction code

Definition at line 200 of file Track3DKalman_module.cc.

{
  rep=0;
  repMC=0;

  // get services
  art::ServiceHandle<geo::Geometry const> geom;
  CLHEP::RandGaussQ gauss(fEngine);

  //////////////////////////////////////////////////////
  // Make a std::unique_ptr<> for the thing you want to put into the event
  // because that handles the memory management for you
  //////////////////////////////////////////////////////
  std::unique_ptr<std::vector<recob::Track> > tcol(new std::vector<recob::Track>);
  std::unique_ptr< std::vector<recob::SpacePoint> >              spcol  (new std::vector<recob::SpacePoint>);
  std::unique_ptr< art::Assns<recob::Track, recob::SpacePoint> > tspassn(new art::Assns<recob::Track, recob::SpacePoint>);
  std::unique_ptr< art::Assns<recob::Track, recob::Cluster> >    tcassn (new art::Assns<recob::Track, recob::Cluster>);
  std::unique_ptr< art::Assns<recob::Track, recob::Hit> >        thassn (new art::Assns<recob::Track, recob::Hit>);

  // define TPC parameters
  TString tpcName = geom->GetLArTPCVolumeName();


  // get input Hit object(s).
  art::Handle< std::vector<recob::Track> > trackListHandle;
  evt.getByLabel(fSpacePtsModuleLabel,trackListHandle);

  art::PtrVector<simb::MCTruth> mclist;

  if (!evt.isRealData())
    {
      art::Handle< std::vector<simb::MCTruth> > mctruthListHandle;
      evt.getByLabel(fGenieGenModuleLabel,mctruthListHandle);

      for (unsigned int ii = 0; ii <  mctruthListHandle->size(); ++ii)
        {
          art::Ptr<simb::MCTruth> mctparticle(mctruthListHandle,ii);
          mclist.push_back(mctparticle);
        }
    }

  //create collection of spacepoints that will be used when creating the Track object
  std::vector< art::Ptr<recob::SpacePoint> > spacepoints;
  art::PtrVector<recob::Track> trackIn;
  // std::cout<<"Run "<<evt.run()<<" Event "<<evt.id().event()<<std::endl;
  mf::LogInfo("Track3DKalman: ") << "There are " <<  trackListHandle->size() << " Track3Dreco/SpacePt tracks/groups (whichever) in this event.";

  art::FindManyP<recob::SpacePoint> fmsp(trackListHandle, evt, fSpacePtsModuleLabel);
  art::FindManyP<recob::Cluster>    fmc (trackListHandle, evt, fSpacePtsModuleLabel);
  art::FindManyP<recob::Hit>        fmh (trackListHandle, evt, fSpacePtsModuleLabel);

  for(unsigned int ii = 0; ii < trackListHandle->size(); ++ii)
    {
      art::Ptr<recob::Track> track(trackListHandle, ii);
      trackIn.push_back(track);

    }

      TVector3 MCOrigin;
      TVector3 MCMomentum;
      // TVector3 posErr(.05,.05,.05); // resolution. 0.5mm
      // TVector3 momErr(.1,.1,0.2);   // GeV
      TVector3 posErr(fPosErr[0],fPosErr[1],fPosErr[2]); // resolution. 0.5mm
      TVector3 momErr(fMomErr[0],fMomErr[1],fMomErr[2]);   // GeV

      // This is strictly for MC
      ///\todo: remove this statement, there is no place for checks on isRealData in reconstruction code
      if (!evt.isRealData())
        {
          // Below breaks are stupid, I realize. But rather than keep all the MC
          // particles I just take the first primary, e.g., muon and only keep its
          // info in the branches of the Ttree. I could generalize later, ...
          for( unsigned int ii = 0; ii < mclist.size(); ++ii )
            {
            //art::Ptr<const simb::MCTruth> mc(mctruthListHandle,i);
            art::Ptr<simb::MCTruth> mc(mclist[ii]);
            for(int jj = 0; jj < mc->NParticles(); ++jj)
              {
                simb::MCParticle part(mc->GetParticle(jj));
                mf::LogInfo("Track3DKalman: ") << "FROM MC TRUTH, the particle's pdg code is: "<<part.PdgCode()<< " with energy = "<<part.E() <<", with energy = "<<part.E()<< " and vtx and momentum in Global (not volTPC) coords are " ;
                MCOrigin.SetXYZ(part.Vx(),part.Vy(),part.Vz()); // V for Vertex
                MCMomentum.SetXYZ(part.Px(),part.Py(),part.Pz());
                MCOrigin.Print();
                MCMomentum.Print();
                repMC = new genf::RKTrackRep(MCOrigin,
                                             MCMomentum,
                                             posErr,
                                             momErr,
                                             part.PdgCode());
                break;
              }
            break;
          }
          //for saving of MC truth
          stMCT->ResizeTo(repMC->getState());
          *stMCT = repMC->getState();
          covMCT-> ResizeTo(repMC->getCov());
          *covMCT = repMC->getCov();
          mf::LogInfo("Track3DKalman: ") <<" repMC, covMC are ... " ;
          repMC->getState().Print();
          repMC->getCov().Print();

        } // !isRealData

      art::PtrVector<recob::Track>::const_iterator trackIter = trackIn.begin();

      nTrks=0;
      while(trackIter!=trackIn.end()) {
        spacepoints.clear();
        spacepoints = fmsp.at(nTrks);

        mf::LogInfo("Track3DKalman: ") << "found "<< spacepoints.size() <<" 3D spacepoint(s).";

        // Add the 3D track to the vector of the reconstructed tracks
        if(spacepoints.size()>0){
          // Insert the GENFIT/Kalman stuff here then make the tracks. Units are cm, GeV.
          const double resolution = 0.5; // dunno, 5 mm
          const int numIT = 3; // 3->1, EC, 6-Jan-2011. Back, 7-Jan-2011.


          //TVector3 mom(0.0,0.0,2.0);
          TVector3 mom(fMomStart[0],fMomStart[1],fMomStart[2]);
          //mom.SetMag(1.);
          TVector3 momM(mom);
          momM.SetX(gauss.fire(momM.X(),momErr.X()/* *momM.X() */));
          momM.SetY(gauss.fire(momM.Y(),momErr.Y()/* *momM.Y() */));
          momM.SetZ(gauss.fire(momM.Z(),momErr.Z()/* *momM.Z() */));
          //std::cout << "Track3DKalman: sort spacepoints by z

          std::sort(spacepoints.begin(), spacepoints.end(), sp_sort_3dz);

          //std::sort(spacepoints.begin(), spacepoints.end(), sp_sort_3dx); // Reverse sort!

          genf::GFFieldManager::getInstance()->init(new genf::GFConstField(0.,0.,0.0));
          genf::GFDetPlane planeG((TVector3)(spacepoints[0]->XYZ()),momM);


          //      std::cout<<"Track3DKalman about to do GAbsTrackRep."<<std::endl;
          // Initialize with 1st spacepoint location and a guess at the momentum.
          rep = new genf::RKTrackRep(//posM-.5/momM.Mag()*momM,
                                     (TVector3)(spacepoints[0]->XYZ()),
                                     momM,
                                     posErr,
                                     momErr,
                                     13);  // mu- hypothesis
          //      std::cout<<"Track3DKalman: about to do GFTrack. repDim is " << rep->getDim() <<std::endl;


          genf::GFTrack fitTrack(rep);//initialized with smeared rep
          // Gonna sort in x cuz I want to essentially transform here to volTPC coords.
          // volTPC coords, cuz that's what the Geant3/Geane stepper wants, as that's its understanding
          // from the Geant4 geometry, which it'll use. EC, 7-Jan-2011.
          int ihit = 0;
          fptsNo = 0;
          for (unsigned int point=0;point<spacepoints.size();++point)
            {

              TVector3 spt3(spacepoints[point]->XYZ());
                if (point > 0 )
                  {
                    double eps (0.1);
                    TVector3 magNew(spt3[0],spt3[1],spt3[2]);
                    TVector3 magLast(spacepoints.back()->XYZ()[0],
                                     spacepoints.back()->XYZ()[1],
                                     spacepoints.back()->XYZ()[2]
                                     );
                    if (!(magNew.Mag()>=magLast.Mag()+eps ||
                          magNew.Mag()<=magLast.Mag()-eps)
                        )
                      continue;
                  }

              if (point%20) // Jump out of loop except on every 20th pt.
                {
                  //continue;
                  // Icarus paper suggests we may wanna decimate our data in order to give
                  // trackfitter a better idea of multiple-scattering. EC, 7-Jan-2011.
                  //if (std::abs(spt3[0]-spacepoints.at(point-1).XYZ()[0]) < 2) continue;
                }
              if (fptsNo<fDimSize)
                {
                  fshx[fptsNo] = spt3[0];
                  fshy[fptsNo] = spt3[1];
                  fshz[fptsNo] = spt3[2];
                }
                fptsNo++;


                MF_LOG_DEBUG("Track3DKalman: ") << "ihit xyz..." << spt3[0]<<","<< spt3[1]<<","<< spt3[2];
                fitTrack.addHit(new genf::PointHit(spt3,resolution),
                                1,//dummy detector id
                                ihit++
                                );
            }

          //      std::cout<<"Track3DKalman about to do GFKalman."<<std::endl;
          genf::GFKalman k;
          //k.setBlowUpFactor(500); // Instead of 500 out of box. EC, 6-Jan-2011.
          //k.setInitialDirection(+1); // Instead of 1 out of box. EC, 6-Jan-2011.
          k.setNumIterations(numIT);
          //      std::cout<<"Track3DKalman back from setNumIterations."<<std::endl;
          try{
            //  std::cout<<"Track3DKalman about to processTrack."<<std::endl;
            k.processTrack(&fitTrack);
            //std::cout<<"Track3DKalman back from processTrack."<<std::endl;
          }
          catch(GFException& e){
            mf::LogError("Track3DKalman: ") << "just caught a GFException."<<std::endl;
            e.what();
            mf::LogError("Track3DKalman: ") << "Exceptions won't be further handled ->exit(1) "<<__LINE__;

            //  exit(1);
          }

          if(rep->getStatusFlag()==0) // 0 is successful completion
            {
              MF_LOG_DEBUG("Track3DKalman: ") << __FILE__ << " " << __LINE__ ;
              MF_LOG_DEBUG("Track3DKalman: ") << "Track3DKalman.cxx: Original plane:";

              if(fGenfPRINT) planeG.Print();
              MF_LOG_DEBUG("Track3DKalman: ") << "Current (fit) reference Plane:";
              if(fGenfPRINT) rep->getReferencePlane().Print();

              MF_LOG_DEBUG("Track3DKalman: ") << "Track3DKalman.cxx: Last reference Plane:";
              if(fGenfPRINT) rep->getLastPlane().Print();

              if(fGenfPRINT)
                {
                  if(planeG!=rep->getReferencePlane())
                    MF_LOG_DEBUG("Track3DKalman: ") <<"Track3DKalman: Original hit plane (not surprisingly) not current reference Plane!"<<std::endl;
                }

              stREC->ResizeTo(rep->getState());
              *stREC = rep->getState();
              covREC->ResizeTo(rep->getCov());
              *covREC = rep->getCov();
              if(fGenfPRINT)
                {
                  MF_LOG_DEBUG("Track3DKalman: ") << " Final State and Cov:";
                  stREC->Print();
                  covREC->Print();
                }
              chi2 = rep->getChiSqu();
              ndf = rep->getNDF();
              nfail = fitTrack.getFailedHits();
              chi2ndf = chi2/ndf;
              TVector3 dircoss = (*trackIter)->VertexDirection<TVector3>();

              for (int ii=0;ii<3;++ii)
                {
                  fpMCT[ii] = MCMomentum[ii]/MCMomentum.Mag();
                  fpREC[ii] = rep->getReferencePlane().getNormal()[ii];
                  fpRECL[ii] = rep->getLastPlane().getNormal()[ii];
                  fpRECt3D[ii] = dircoss[ii];
                }
              fpMCT[3] = MCMomentum.Mag();
              fpREC[3] = -1.0/(*stREC)[0][0];

              evtt = (unsigned int) evt.id().event();
              mf::LogInfo("Track3DKalman: ") << "Track3DKalman about to do tree->Fill(). Chi2/ndf is " << chi2/ndf << ". All in volTPC coords .... pMCT[0-3] is " << fpMCT[0] << ", " << fpMCT[1] << ", " << fpMCT[2] << ", " << fpMCT[3] << ". pREC[0-3] is " << fpREC[0] << ", "<< fpREC[1] << ", " << fpREC[2] << ", " << fpREC[3] << ".";

              tree->Fill();

              // Get the clusters associated to each track in induction and collection view
              std::vector< art::Ptr<recob::Cluster> > clusters = fmc.at(nTrks);
              std::vector< art::Ptr<recob::Hit> > hits = fmh.at(nTrks);

              // Use new Track constructor... EC, 21-Apr-2011.
              std::vector<TVector3> dircos(spacepoints.size());
              dircos[0] = TVector3(fpREC);
              dircos.back() = TVector3(fpRECL);

              // fill a vector of TVector3 with the space points used to make this track
              std::vector<TVector3> xyz(spacepoints.size());
              size_t spStart = spcol->size();
              for(size_t tv = 0; tv < spacepoints.size(); ++tv){
                xyz[tv] = TVector3(spacepoints[tv]->XYZ());
                spcol->push_back(*(spacepoints[tv].get()));
              }
              size_t spEnd = spcol->size();

              tcol->push_back(recob::Track(recob::TrackTrajectory(recob::tracking::convertCollToPoint(xyz),
                                                                  recob::tracking::convertCollToVector(dircos),
                                                                  recob::Track::Flags_t(xyz.size()), false),
                                           0, -1., 0, recob::tracking::SMatrixSym55(), recob::tracking::SMatrixSym55(), tcol->size()-1));

              // associate the track with its clusters and tracks
              util::CreateAssn(*this, evt, *tcol, clusters, *tcassn);
              util::CreateAssn(*this, evt, *tcol, hits,     *thassn);

              // associate the track to the space points
              util::CreateAssn(*this, evt, *tcol, *spcol, *tspassn, spStart, spEnd);

            } // getStatusFlag
        } // spacepoints.size()>0

        //
        //std::cout<<"Track3DKalman found "<< tcol->size() <<" 3D track(s)"<<std::endl;
        if(trackIter!=trackIn.end()) trackIter++;

        nTrks++;

      } // end loop over Track3Dreco/SpacePt tracks/groups (whichever) we brought into this event.

      evt.put(std::move(tcol));
      evt.put(std::move(spcol));
      evt.put(std::move(tcassn));
      evt.put(std::move(thassn));
      evt.put(std::move(tspassn));

} // end method

Member Data Documentation

Float_t trkf::Track3DKalman::chi2

private

Definition at line 82 of file Track3DKalman_module.cc.

Float_t trkf::Track3DKalman::chi2ndf

private

Definition at line 83 of file Track3DKalman_module.cc.

TMatrixT<Double_t>* trkf::Track3DKalman::covMCT

private

Definition at line 79 of file Track3DKalman_module.cc.

TMatrixT<Double_t>* trkf::Track3DKalman::covREC

private

Definition at line 81 of file Track3DKalman_module.cc.

unsigned int trkf::Track3DKalman::evtt

private

Definition at line 91 of file Track3DKalman_module.cc.

unsigned int trkf::Track3DKalman::fDimSize

private

Definition at line 97 of file Track3DKalman_module.cc.

CLHEP::HepRandomEngine& trkf::Track3DKalman::fEngine

private

Definition at line 104 of file Track3DKalman_module.cc.

std::string trkf::Track3DKalman::fG4ModuleLabel

private

Definition at line 73 of file Track3DKalman_module.cc.

bool trkf::Track3DKalman::fGenfPRINT

private

Definition at line 74 of file Track3DKalman_module.cc.

std::string trkf::Track3DKalman::fGenieGenModuleLabel

private

Definition at line 72 of file Track3DKalman_module.cc.

std::vector<double> trkf::Track3DKalman::fMomErr

private

Definition at line 100 of file Track3DKalman_module.cc.

std::vector<double> trkf::Track3DKalman::fMomStart

private

Definition at line 101 of file Track3DKalman_module.cc.

Float_t* trkf::Track3DKalman::fpMCT

private

Definition at line 88 of file Track3DKalman_module.cc.

std::vector<double> trkf::Track3DKalman::fPosErr

private

Definition at line 99 of file Track3DKalman_module.cc.

Float_t* trkf::Track3DKalman::fpREC

private

Definition at line 87 of file Track3DKalman_module.cc.

Float_t* trkf::Track3DKalman::fpRECL

private

Definition at line 86 of file Track3DKalman_module.cc.

Float_t* trkf::Track3DKalman::fpRECt3D

private

Definition at line 85 of file Track3DKalman_module.cc.

unsigned int trkf::Track3DKalman::fptsNo

private

Definition at line 93 of file Track3DKalman_module.cc.

Float_t* trkf::Track3DKalman::fshx

private

Definition at line 94 of file Track3DKalman_module.cc.

Float_t* trkf::Track3DKalman::fshy

private

Definition at line 95 of file Track3DKalman_module.cc.

Float_t* trkf::Track3DKalman::fshz

private

Definition at line 96 of file Track3DKalman_module.cc.

std::string trkf::Track3DKalman::fSpacePtsModuleLabel

private

Definition at line 71 of file Track3DKalman_module.cc.

int trkf::Track3DKalman::ndf

private

Definition at line 90 of file Track3DKalman_module.cc.

int trkf::Track3DKalman::nfail

private

Definition at line 89 of file Track3DKalman_module.cc.

unsigned int trkf::Track3DKalman::nTrks

private

Definition at line 92 of file Track3DKalman_module.cc.

genf::GFAbsTrackRep* trkf::Track3DKalman::rep

private

Definition at line 103 of file Track3DKalman_module.cc.

genf::GFAbsTrackRep* trkf::Track3DKalman::repMC

private

Definition at line 102 of file Track3DKalman_module.cc.

TMatrixT<Double_t>* trkf::Track3DKalman::stMCT

private

Definition at line 78 of file Track3DKalman_module.cc.

TMatrixT<Double_t>* trkf::Track3DKalman::stREC

private

Definition at line 80 of file Track3DKalman_module.cc.

TTree* trkf::Track3DKalman::tree

private

Definition at line 76 of file Track3DKalman_module.cc.

---

The documentation for this class was generated from the following file:

• Track3DKalman_module.cc