TrackInfo_module.cc
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////
2 /// \brief Track information
3 /// \author Leo Aliaga laliaga@fnal.gov
4 ////////////////////////////////////////////////////////
5 #include <string>
6 
11 #include "art_root_io/TFileService.h"
12 #include "fhiclcpp/ParameterSet.h"
14 
15 #include "TrackInfo/TrackInfoFxs.h"
16 #include "TrackInfo/TrackInfoObj.h"
17 #include "Geometry/Geometry.h"
18 #include "Geometry/LiveGeometry.h"
19 #include "MCCheater/BackTracker.h"
20 #include "NovaDAQConventions/DAQConventions.h"
21 #include "RecoBase/CellHit.h"
22 #include "RecoBase/Cluster.h"
23 #include "RecoBase/Track.h"
24 #include "TimingFit/TimingFitAlg.h"
25 #include "Utilities/AssociationUtil.h"
27 #include "CVN/func/Result.h"
30 
31 #include "TFile.h"
32 #include "TMVA/Reader.h"
33 
34 namespace trackinfo {
35 
36  class TrackInfo : public art::EDProducer {
37  public:
38  explicit TrackInfo(fhicl::ParameterSet const & pset);
39  virtual ~TrackInfo();
40 
41  void produce (art::Event & evt);
42 
43  protected:
44  void Init();
45 
47  bpfit::dEdxCalculator fdEdxCalc; ///< helper class for computing dEdx values
48  double fdxTOL; ///< cut off value for dx, if dx is less than this, don't compute dE/dx
52  };
53 }
54 
55 
56 namespace trackinfo
57 {
58 
59  //----------------------------------------------------------------------
61  : EDProducer(pset),
63  fdxTOL (pset.get<double> ("dxTOL")),
64  fSlicerLabel (pset.get< std::string >("SlicerLabel")),
65  fTrackLabel (pset.get< std::string >("TrackLabel")){
66  produces< std::vector<trackinfo::TrackInfoObj> >();
67  produces< art::Assns<trackinfo::TrackInfoObj, rb::Track> >();
68  }
69 
70 
71  //----------------------------------------------------------------------
73  if(fTrackInfoFxs) delete fTrackInfoFxs;
74  }
75 
76  //-------------------------------------------------------------------
77  void TrackInfo::Init(){
79  }
80 
81  //-------------------------------------------------------------------
83 
84  Init();
89  evt.getByLabel(fSlicerLabel,slicevec);
91 
92  if(slicevec->empty()) {
93  mf::LogWarning ("No Slices")<<"No Slices in the input file";
94  return;
95  }
96  art::FindManyP<rb::Track> trackAssnList(slicevec, evt, fTrackLabel);
97  std::unique_ptr< std::vector<trackinfo::TrackInfoObj> > outputObjects(new std::vector<trackinfo::TrackInfoObj>);
98  std::unique_ptr< art::Assns<trackinfo::TrackInfoObj, rb::Track> > assoc(new art::Assns<trackinfo::TrackInfoObj, rb::Track>);
99 
100  for(size_t i = 0; i < slicevec->size(); ++i){
101  art::Ptr<rb::Cluster> slice(slicevec,i);
102 
103  if(slice->IsNoise()) {
104  continue;
105  }
106 
107  const std::vector< art::Ptr<rb::Track> > sliceTracks = trackAssnList.at(i);
108  for(size_t iTrack = 0; iTrack < sliceTracks.size(); ++iTrack){
109  trackinfo::TrackInfoObj trackinfoobj;
110  art::Ptr<rb::Track> track = sliceTracks[iTrack];
111 
112  // ignore bad track
113  double tdx = track->Dir().X();
114  double tdy = track->Dir().Y();
115  double tdz = track->Dir().Z();
116  bool bad_track_dir = (tdx < -1) || (tdx > 1) ||
117  (tdy < -1) || (tdy > 1) ||
118  (tdz < -1) || (tdz > 1);
119  if(bad_track_dir){
120  outputObjects->push_back(trackinfoobj);
121  util::CreateAssn(evt,*(outputObjects.get()),track,*(assoc.get()));
122  continue;
123  }
124 
125  //We are only interested on the dE,dx and s info from dEdxCalculator:
126  const int fake_pdg = 13;
127  fdEdxCalc.computeDEDX(track,fake_pdg);
128  std::vector<double> dE, dx, s;
129  fdEdxCalc.getDE(dE,fdxTOL);
130  fdEdxCalc.getDX(dx,fdxTOL);
131  fdEdxCalc.getS( s ,fdxTOL);
132  trackinfoobj.SetAvedEdxTrackLast10cm( fTrackInfoFxs->getAveTrackdEdxEnd(dE, dx, s, track->TotalLength(), 10.) );
133  trackinfoobj.SetAvedEdxTrackLast20cm( fTrackInfoFxs->getAveTrackdEdxEnd(dE, dx, s, track->TotalLength(), 20.) );
134  trackinfoobj.SetAvedEdxTrackLast30cm( fTrackInfoFxs->getAveTrackdEdxEnd(dE, dx, s, track->TotalLength(), 30.) );
135  trackinfoobj.SetAvedEdxTrackLast40cm( fTrackInfoFxs->getAveTrackdEdxEnd(dE, dx, s, track->TotalLength(), 40.) );
136  outputObjects->push_back(trackinfoobj);
137  util::CreateAssn(evt,*(outputObjects.get()),track,*(assoc.get()));
138 
139  } // end of loop over tracks on slice
140 
141  } //loop over slices
142 
143  evt.put(std::move(outputObjects));
144  evt.put(std::move(assoc));
145 
146  } //end produce
147 
148 } //end namespace
149 
back track the reconstruction to the simulation
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
void SetAvedEdxTrackLast20cm(float AvedEdxTrackLast20cm_val)
Definition: TrackInfoObj.h:25
void SetAvedEdxTrackLast30cm(float AvedEdxTrackLast30cm_val)
Definition: TrackInfoObj.h:26
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
Definition: event.h:19
bpfit::dEdxCalculator fdEdxCalc
helper class for computing dEdx values
void getS(std::vector< double > &stemp, double dxTol)
void computeDEDX(art::Ptr< rb::Track > &track, int pdg=13)
Compute dE/dx for all cells along this track and the fitsum object that went with it...
double fdxTOL
cut off value for dx, if dx is less than this, don&#39;t compute dE/dx
DEFINE_ART_MODULE(TestTMapFile)
double dE
void SetAvedEdxTrackLast40cm(float AvedEdxTrackLast40cm_val)
Definition: TrackInfoObj.h:27
Calculates dE/dx in cells passed through by a track.
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
const XML_Char * s
Definition: expat.h:262
trackinfo::TrackInfoFxs * fTrackInfoFxs
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
double dx[NP][NC]
Result for CVN.
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
int evt
TrackInfo(fhicl::ParameterSet const &pset)
float getAveTrackdEdxEnd(std::vector< double > &dEtemp, std::vector< double > &dxtemp, std::vector< double > &stemp, double trk_len, double wanted_dist)
void SetAvedEdxTrackLast10cm(float AvedEdxTrackLast10cm_val)
Definition: TrackInfoObj.h:24
Track information.
Definition: FillReco.h:15
void produce(art::Event &evt)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
std::string fCosmicTrackLabel
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
void getDE(std::vector< double > &dEtemp, double dxTol)
...can also get the individual results...
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
Encapsulate the geometry of one entire detector (near, far, ndos)
void getDX(std::vector< double > &dxtemp, double dxTol)
enum BeamMode string