QeFinder_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////
2 /// \brief A module for finding numu CC QE interactions
3 /// \author Nicholas Raddatz - raddatz@physics.umn.edu
4 /// \date February 2013
5 ////////////////////////////////////////////////////////////////////
6 
7 //ROOT includes
8 #include "TMVA/Reader.h"
9 
10 // Framework includes
13 #include "fhiclcpp/ParameterSet.h"
15 
16 //NOvA includes
17 #include "QEEventFinder/QePId.h"
18 #include "ReMId/ReMId.h"
19 #include "RecoBase/Track.h"
20 #include "NumuEnergy/NumuE.h"
21 #include "Utilities/AssociationUtil.h"
23 #include "Geometry/Geometry.h"
24 
25 
26 namespace qeef{
27  /// A QE Event Finder for muons
28  class QeFinder : public art::EDProducer{
29  public:
30  explicit QeFinder(fhicl::ParameterSet const& pset);
31  virtual ~QeFinder();
32 
33  void beginJob();
34  void produce(art::Event& evt);
35 
36  private:
37  void Init();
38 
47  TMVA::Reader *knn1trk, *knn2trk;
48  Float_t offEK, ediff, ediffZ, dedx;
49 
50  };
51 }
52 
53 namespace qeef{
54 
55  //.......................................................................
57  knn1trk(nullptr),
58  knn2trk(nullptr)
59  {
60  fSliceLabel = (pset.get< std::string >("SliceLabel"));
61  fTrackLabel = (pset.get< std::string >("TrackLabel"));
62  fPIDLabel = (pset.get< std::string >("PIDLabel"));
63  fEnergyLabel = (pset.get< std::string >("EnergyLabel"));
64  fKnnLibPath = (pset.get< std::string >("KnnLibPath"));
65  produces< std::vector<qeef::QePId> >();
66  produces< art::Assns<qeef::QePId, rb::Cluster> >();
67  }
68 
69  //.......................................................................
71  {
72  if (knn1trk) delete knn1trk;
73  if (knn2trk) delete knn2trk;
74  }
75 
76  //.......................................................................
78  {
79  trainfileloaded = false;
80  }
81 
82  //.......................................................................
84  {
85  if(trainfileloaded) return;
86 
88  // Expand fcl input lib path
90  // set up FD training file
91  if(geom->DetId() == novadaq::cnv::kFARDET){
92  fknn1trkfile = knnlibpath+"TMVAClassification_KNN_qe_1track.weights.xml";
93  fknn2trkfile = knnlibpath+"TMVAClassification_KNN_qe_2track.weights.xml";
94  }
95  // set up ND training file
96  else if(geom->DetId() == novadaq::cnv::kNEARDET){
97  fknn1trkfile = knnlibpath+"TMVAClassification_KNN_qe_1track.nd.weights.xml";
98  fknn2trkfile = knnlibpath+"TMVAClassification_KNN_qe_2track.nd.weights.xml";
99  }
100  else{ std::cerr << "This file isn't near or far detector, algorithm doesn't work on any other" << std::endl;}
101 
102  knn1trk = new TMVA::Reader();
103  knn1trk->AddVariable("offEK",&offEK);
104  knn1trk->AddVariable("ediff",&ediff);
105  knn1trk->AddVariable("ediffZ",&ediffZ);
106 
107  knn1trk->BookMVA("KNN method",fknn1trkfile);
108 
109  knn2trk = new TMVA::Reader();
110  knn2trk->AddVariable("offEK",&offEK);
111  knn2trk->AddVariable("ediff",&ediff);
112  knn2trk->AddVariable("ediffZ",&ediffZ);
113  knn2trk->AddVariable("dedx",&dedx);
114 
115  knn2trk->BookMVA("KNN method",fknn2trkfile);
116 
117  trainfileloaded = true;
118  }
119 
120  //.......................................................................
122  {
123  Init();
124 
125  // going to produce qe pid objects and associations with slices
126  std::unique_ptr< std::vector<qeef::QePId> > qeefOut(new std::vector<qeef::QePId>);
127  std::unique_ptr< art::Assns<qeef::QePId, rb::Cluster> > qeAss(new art::Assns<qeef::QePId, rb::Cluster>);
128 
129 
130  // get all of slices out of the event
132  evt.getByLabel(fSliceLabel,slices);
133 
134  // set up associations between slices and tracks
135  art::FindManyP<rb::Track> trkAss(slices,evt,fTrackLabel);
136  // set up associations between slices and energy
138 
139  for(size_t isl =0; isl < slices->size(); ++isl){
140 
141  // Get the slice
142  art::Ptr<rb::Cluster> slice(slices,isl);
143  // skip if it is noise
144  if(slice->IsNoise() || !trkAss.isValid()){ continue; }
145 
146  // make a qepid object with default pid value of -5
147  qeef::QePId qepid(14,-5);
148 
149  // get the tracks that are associated with this slice
150  const std::vector<art::Ptr<rb::Track> > tracks = trkAss.at(isl);
151 
152  // if more than 3 or no tracks no qe pid will be made so skip this
153  if(tracks.size() < 3 && tracks.size() > 0){
154 
155  // set up associations between tracks and remid
156  art::FindOneP<remid::ReMId> remAss(tracks,evt,fPIDLabel);
157 
158  //Find index of highest pid track
159  unsigned int bestId = remid::HighestPIDTrack(tracks, fPIDLabel, evt);
160 
161  int n3dtrk = 0; // keep track of the number of 3d tracks
162  std::vector<art::Ptr<rb::CellHit> > trkHits; // cellhits from every reco track
163 
164  // loop over tracks and find the highest pid one
165  for(size_t itrk = 0; itrk < tracks.size(); ++itrk){
166  // get the track
167  const art::Ptr<rb::Track> track = tracks[itrk];
168 
169  if(track->Is3D()){ ++n3dtrk; }
170 
171  // get all of the cell hits associated with this track
172  for(unsigned int ihit = 0; ihit < track->NCell(); ++ihit){
173  trkHits.push_back(track->Cell(ihit));
174  }
175 
176  }
177 
178  // continue if highest pid track found
179  if(bestId != 999){
180 
181  // check number of 3d trks
182  if((n3dtrk == 1 && tracks.size() == 1) ||
183  (n3dtrk >=1 && tracks.size() == 2)){
184 
185  // get all the hits not assocaited with a track in the event
186  rb::Cluster nontrkHits;
187  for(unsigned int slhit = 0; slhit < slice->NCell(); ++slhit){
188  bool intrk = false;
189  // check if this hit is in any track
190  for(unsigned int trhit = 0; trhit < trkHits.size(); ++trhit){
191  if(slice->Cell(slhit) == trkHits[trhit]){
192  intrk = true;
193  break;
194  }
195  }
196  if(!intrk){ nontrkHits.Add(slice->Cell(slhit)); }
197  }
198 
199  // knn variables
200  float offE;
201  float ediff;
202  float ediffZ;
203  float dedx;
204 
205  // get muon pid and track
206  art::Ptr<remid::ReMId> murem = remAss.at(bestId);
207  art::Ptr<rb::Track> mutrk = tracks[bestId];
208 
209  // get numu energy estimates
210  if(eAss.isValid()){
211  std::vector<art::Ptr<numue::NumuE> > numuEnergy = eAss.at(isl);
212  if(numuEnergy.size() > 0){
213 
214  // get the energy variables for the qe pid
215  float angleQEE = numuEnergy[0]->AngleQEE();
216  float errAngleQEE = numuEnergy[0]->AngleQEError();
217  float trkQEE = numuEnergy[0]->TrkQEE();
218 
219  // Don't pid for cases where angleQEE or trkqeE is nonsense
220  if( angleQEE < 0.0 || angleQEE > 40 || errAngleQEE < 0.0 || errAngleQEE > 40 || trkQEE > 20){
221  qepid.SetNtrk(tracks.size());
222  qepid.SetVal(-1);
223  }
224  else{
225 
226  ediff = (angleQEE - trkQEE)/trkQEE;
227  ediffZ = (angleQEE - trkQEE)/errAngleQEE;
228 
229  if(tracks.size() == 1){
230  double vertE = 0.0;
231  if(nontrkHits.NCell() > 0){
232  vertE = nontrkHits.TotalGeV();
233  }
234  offE = vertE/mutrk->TotalGeV();
235 
236  // gather all pid variables
237  std::vector<float> knnInputs(3);
238  knnInputs[0] = offE;
239  knnInputs[1] = ediff;
240  knnInputs[2] = ediffZ;
241  double kNN = knn1trk->EvaluateMVA(knnInputs,"KNN method");
242 
243  // update qe pid info
244  qepid.SetNtrk(1);
245  qepid.SetOffTrkE(offE);
246  qepid.SetEdiff(ediff);
247  qepid.SetEdiffZ(ediffZ);
248  qepid.SetVal(kNN);
249  } // end of 1 track case
250  else{
251  int nonmuId = 1;
252  if(bestId == 1){ nonmuId = 0; }
253  art::Ptr<remid::ReMId> nonmurem = remAss.at(nonmuId);
254  art::Ptr<rb::Track> nonmutrk = tracks[nonmuId];
255  double vertE = 0.0;
256  if(nontrkHits.NCell() > 0){
257  vertE = nontrkHits.TotalGeV();
258  }
259  offE = vertE/(mutrk->TotalGeV()+nonmutrk->TotalGeV());
260  dedx = nonmurem->AvgDedx()/murem->AvgDedx();
261 
262  // gather all pid variables
263  std::vector<float> knnInputs(4);
264  knnInputs[0] = offE;
265  knnInputs[1] = ediff;
266  knnInputs[2] = ediffZ;
267  knnInputs[3] = dedx;
268  double kNN = knn2trk->EvaluateMVA(knnInputs,"KNN method");
269 
270  // update qe pid info
271  qepid.SetNtrk(2);
272  qepid.SetOffTrkE(offE);
273  qepid.SetEdiff(ediff);
274  qepid.SetEdiffZ(ediffZ);
275  qepid.SetDedx(dedx);
276  qepid.SetVal(kNN);
277  } // end of 2 track case
278 
279  } // end of else nonsense variables
280  } // check if we have numu energy for this slice
281  } // check that numu energy association is valid
282 
283  } // check if we have the correct number of tracks for qe piding
284 
285  } // check if we found the best muon
286 
287 
288 
289  } // end of if < 3 total tracks
290 
291  // add this pid object to the vector of ones to add to the event
292  qeefOut->push_back(qepid);
293  // make an association to the slice it came from
294  util::CreateAssn(*this,evt,*(qeefOut.get()),slice,*(qeAss.get()));
295 
296  } // end loop over slices
297 
298  evt.put(std::move(qeefOut));
299  evt.put(std::move(qeAss));
300 
301  } // end of produce
302 
303 } // end of qeef namespace
304 
305 namespace qeef
306 {
308 }
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.
std::string fknn1trkfile
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
TMVA::Reader * knn2trk
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
void SetEdiff(double diff)
Definition: QePId.cxx:51
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
virtual ~QeFinder()
OStream cerr
Definition: OStream.cxx:7
void produce(art::Event &evt)
std::string fTrackLabel
DEFINE_ART_MODULE(TestTMapFile)
std::string fknn2trkfile
std::string fPIDLabel
double AvgDedx() const
Definition: ReMId.cxx:68
std::string fKnnLibPath
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
QeFinder(fhicl::ParameterSet const &pset)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void SetEdiffZ(double diffZ)
Definition: QePId.cxx:57
Far Detector at Ash River, MN.
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
void SetDedx(double Dedx)
Definition: QePId.cxx:63
void SetOffTrkE(double E)
Definition: QePId.cxx:45
Near Detector in the NuMI cavern.
A QE Event Finder for muons.
TMVA::Reader * knn1trk
void SetNtrk(int ntrks)
Definition: QePId.cxx:39
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void geom(int which=0)
Definition: geom.C:163
std::string fSliceLabel
A module for finding numu CC QE interactions.
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
std::string fEnergyLabel
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
void SetVal(double val)
Definition: PID.h:24
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
Definition: Prong.h:71
unsigned int HighestPIDTrack(const std::vector< art::Ptr< rb::Track > > &sliceTracks, const std::string &remidModuleLabel, const art::Event &e)
Definition: ReMId.cxx:249