TrackEva_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    TrackEva.cxx
//
// \brief   Generating a root file containing information you need for getting 
//          track reconstruction efficiency.
//
// \version $Id: TrackEva.cxx,v 1.3 2012-09-25 06:30:12 gsdavies Exp $
// \author  ricken@fnal.gov
////////////////////////////////////////////////////////////////////////
#include "RecoBase/Track.h"
#include "Simulation/Particle.h"
#include "Simulation/ParticleNavigator.h"
#include "nusimdata/SimulationBase/MCTrajectory.h"
#include "MCCheater/BackTracker.h"


// Framework includes
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/TFileService.h"

#include "TVector3.h"
#include "TLorentzVector.h"
#include "TNtuple.h"
#include "TH1I.h"

#include <string>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <fstream>


class TH1I;
class TNtuple;

namespace simb { class MCTrajectory; }

namespace vf {
  class TrackEva : public art::EDAnalyzer
    {
    public:
      explicit TrackEva(fhicl::ParameterSet const &pset);
      virtual ~TrackEva();
                        
      void analyze(const art::Event &evt);
      void beginJob();
      void reconfigure(fhicl::ParameterSet const &pset);

      static const int GAMMAID;
      static const int PIZEROID;
      static const int PIPLUSID;
      static const int PIMINUSID;
      static const int NUEID;
      static const int NUMUID;
      static const int NUTAUID;
      static const int MUONID;
      static const int PROTONID;

      
    private:

      double Tracklength(const simb::MCTrajectory & trajectory);

      TNtuple*    PROTON_RC;
      TNtuple*    MUON_RC;
      TNtuple*    PROTON_MC;
      TNtuple*    MUON_MC;
      
      TH1I*       track_Nb;
      TH1I*       proton_Nb;
      TH1I*       proton_Tr;
      TH1I*       muon_Nb;
      TH1I*       muon_Tr;

      std::string Track_RC;
      double      fMatching;
      double      fEdgeX;
      double      fEdgeY;
      double      fEdgeZ;
      double      fVTX;
      double      fVTY;
      double      fVTZ;
    };


  const int TrackEva::GAMMAID = 22;
  const int TrackEva::PIZEROID = 111;
  const int TrackEva::PIPLUSID = 211;
  const int TrackEva::PIMINUSID = -211;
  const int TrackEva::NUEID = 12;
  const int TrackEva::MUONID = 13;
  const int TrackEva::NUMUID = 14;
  const int TrackEva::NUTAUID = 16;
  const int TrackEva::PROTONID = 2212;


  TrackEva::TrackEva(fhicl::ParameterSet const& pset) 
    : EDAnalyzer(pset)
  {
      reconfigure(pset);
  }


//......................................................................
  void TrackEva::reconfigure(fhicl::ParameterSet const &pset)
  {
    Track_RC        =pset.get< std::string >("Track_RC" );

    fMatching       =pset.get< double      >("MatchingPurity" );

    fEdgeX          =pset.get< double      >("EdgeX");
    fEdgeY          =pset.get< double      >("EdgeY");
    fEdgeZ          =pset.get< double      >("EdgeZ");

  }

//......................................................................

  TrackEva::~TrackEva(){}
  

//......................................................................

  void TrackEva::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;
   

    //The following two ntuples can give info: 
    // efficiency VS track length 
 
    MUON_MC = tfs->make<TNtuple>("MUON_MC","muon_mc_info",
                                   "energy:trackLength:Dir_Cos:ID");
   
    PROTON_MC = tfs->make<TNtuple>("PROTON_MC","proton_mc_info",
                                   "energy:trackLength:Dir_Cos:ID");

    MUON_RC = tfs->make<TNtuple>("MUON_RC","muon_rc_info",
                                   "tracklength_diff:energy_diff:Angle:TL_MC:E_MC:Dir_Cos:ID");
   
    PROTON_RC = tfs->make<TNtuple>("PROTON_RC","proton_rc_info",
                                   "tracklength_diff:energy_diff:Angle:TL_MC:E_MC:Dir_Cos:ID");
 

    track_Nb = tfs->make<TH1I>("track_Nb","Total Reconstructed Track-Number",
                                          1001,-0.5,1000.5); 

    proton_Nb =  tfs->make<TH1I>("proton_Nb","Total Proton Number",
                                             501,-0.5,500.5);
    proton_Tr =  tfs->make<TH1I>("proton_Tr","Reco Proton Track Number",
                                             500,-0.5,500.5);
    
    muon_Nb =   tfs->make<TH1I>("muon_Nb","Total Muon Number",
                                           1001,-0.5,1000.5);
    muon_Tr =   tfs->make<TH1I>("muon_Tr","Reco Muon Track Number",
                                           1001,-0.5,1000.5);

  }

  void TrackEva::analyze(const art::Event& evt)
  {
    std::vector<double>  Mu_E, P_E, LMC_mu, LMC_p;      
    std::vector<TVector3> P_Momenta, Mu_Momenta;
  
    TVector3 momentum_P, momentum_Mu;
    
    art::Handle< std::vector<rb::Track> > trackhandle;
    try{ evt.getByLabel(Track_RC,trackhandle); }
    catch(...){ std::cout<<"couldn't get the trackhandle "<<std::endl; return;}

    //make a backtracker
    art::ServiceHandle<cheat::BackTracker> bt;
    const sim::ParticleNavigator& pnav = bt->ParticleNavigator();

 
    //==============================================================

    int NbTrack = trackhandle->size();
    track_Nb->Fill(NbTrack);

    // Get MC infor of the proton and muon:

    std::vector<int> MuonID, ProtonID;
 
    for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i){
      const sim::Particle *p = (*i).second;

      if (p->PdgCode()==PROTONID ) {
        LMC_p.push_back( Tracklength(p->Trajectory()));
        momentum_P.SetXYZ(p->Px(),
                          p->Py(),
                          p->Pz());
        P_Momenta.push_back(momentum_P);
        P_E.push_back(p->E());
        ProtonID.push_back(p->TrackId());
      }
      else if (p->PdgCode()==MUONID ){
        LMC_mu.push_back( Tracklength(p->Trajectory())); 
        momentum_Mu.SetXYZ(p->Px(),
                           p->Py(),
                           p->Pz());
        Mu_Momenta.push_back(momentum_Mu);
        Mu_E.push_back( p->E());
        MuonID.push_back(p->TrackId());
      }
    }

    //Filling proton mc ntuple
    for (unsigned int n=0; n!=P_E.size();n++){
      PROTON_MC->Fill(P_E[n],
                      LMC_p[n],
                      P_Momenta[n].Z()/P_Momenta[n].Mag(),
                      ProtonID[n]);
    }

    //Filling muon mc ntuple
    for (unsigned int n=0;n!=Mu_E.size();++n){
      MUON_MC->Fill(Mu_E[n],
                    LMC_mu[n],
                    Mu_Momenta[n].Z()/Mu_Momenta[n].Mag(),
                    MuonID[n]);
    }

    //Filling MC truth particle number info:

    proton_Nb->Fill(P_E.size());
    muon_Nb  ->Fill(Mu_E.size());


    // Let's count qualified Tracks:

    int mu_Track=0;
    int p_Track=0;
  
    for (int i=0; i!=NbTrack; i++){
      art::Ptr<rb::Track> track(trackhandle, i);
      int cell_Nb = track->NCell();
      
      //Muon Track?
      for (unsigned int j=0; j!=MuonID.size(); j++){
        int cell_count=0;

        for (int idx = 0; idx != cell_Nb; idx++){

          bool real_hit= false;

          std::vector<cheat::TrackIDE> trackides= bt->HitToTrackIDE(track->Cell(idx));

          for (unsigned int pts=0; pts!=trackides.size(); pts++){
            if (trackides[pts].trackID==MuonID[j]) real_hit=true;
          }

          if (real_hit) cell_count++;
        }
        if ((cell_count+0.0)/(cell_Nb+0.0)>=fMatching){
          mu_Track++;
          MUON_RC->Fill(
                        track->TotalLength()-LMC_mu[j],
                        track->TotalGeV()-Mu_E[j],
                        Mu_Momenta[j].Angle(track->Dir()),
                        LMC_mu[j],
                        Mu_E[j],
                        Mu_Momenta[j].Z()/Mu_Momenta[j].Mag(),
                        MuonID[j]
                        );               
        }
      }     


      //Proton Track?
      for (unsigned int j=0; j<ProtonID.size(); j++){
        int cell_count=0;
       
        for (int idx = 0; idx < cell_Nb; idx++){

          bool real_hit= false;
          std::vector<cheat::TrackIDE> trackides= bt->HitToTrackIDE(track->Cell(idx));
          for (unsigned int pts=0; pts!=trackides.size(); pts++){
            if (trackides[pts].trackID==ProtonID[j]) real_hit=true;
          }

          if (real_hit) cell_count++;
        }
        if ((cell_count+0.0)/(cell_Nb+0.0)>fMatching){
          p_Track++;
          PROTON_RC->Fill(
                          track->TotalLength()-LMC_p[j],
                          track->TotalGeV()-P_E[j],
                          P_Momenta[j].Angle(track->Dir()),
                          LMC_p[j],
                          P_E[j],
                          P_Momenta[j].Z()/P_Momenta[j].Mag(),
                          ProtonID[j]
                          );               
        }
      }
    }


    proton_Tr -> Fill(p_Track);
    muon_Tr   -> Fill(mu_Track);
    
    return;
  }

  //Calculate the track length inside the detector
  double TrackEva::Tracklength(const simb::MCTrajectory& trajectory){

    double length=0;

    for(simb::MCTrajectory::const_iterator it = trajectory.begin(); it!=trajectory.end()-1; ++it)
      {
        TLorentzVector pos1 = it->first;
        TLorentzVector pos2 = (it+1)->first;
          
        if (fabs(pos2.X())<fEdgeX && fabs(pos2.Y())<fEdgeY){
          if(pos2.Z()>9.96604 && pos2.Z()<fEdgeZ){
            length+=(pos1.Vect()-pos2.Vect()).Mag();
          }
        }
        
      }
    return (length);
    
  }


  DEFINE_ART_MODULE(TrackEva)
   
}