SPIDBuilder_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Build and write out the ShowerPID objects and make shower PID
3 ///
4 /// \todo Jianming's oscillation weighting: Need to find standard
5 // infrastructure to do this and replace his
6 /// code, which I stripped out of this class
7 ///
8 /// \author Alex Smith smith@physics.umn.edu
9 ////////////////////////////////////////////////////////////////////////
10 
11 // ROOT includes
12 #include "TFile.h"
13 #include "TH1D.h"
14 #include "TLorentzVector.h"
15 #include "TMath.h"
16 #include "TMultiLayerPerceptron.h"
17 #include "TRandom.h"
18 #include "TTree.h"
19 
20 // NOvA includes
21 #include "Geometry/Geometry.h"
22 #include "RecoBase/Cluster.h"
23 #include "RecoBase/PID.h"
24 #include "RecoBase/FilterList.h"
25 #include "RecoBase/Vertex.h"
26 #include "Utilities/AssociationUtil.h"
27 #include "Calibrator/Calibrator.h"
28 
29 // ART includes
33 #include "art_root_io/TFileService.h"
41 
42 
43 #include "fhiclcpp/ParameterSet.h"
46 #include "ShowerLID/NuEEnergyAlg.h"
47 #include "ShowerLID/SPIDAlg.h"
48 #include "ShowerLID/ShowerLID.h"
50 
51 #include <cmath>
52 
53 namespace slid {
54 
55  class SPIDBuilder : public art::EDProducer {
56  public:
57  explicit SPIDBuilder(const fhicl::ParameterSet& pset);
58  ~SPIDBuilder();
59 
60  virtual void produce(art::Event& evt);
61 
62  virtual void reconfigure(const fhicl::ParameterSet& pset);
63 
64  //virtual void beginJob();
65  virtual void beginRun(art::Run& run);
66  //virtual void endJob();
67 
68  protected:
72  std::vector<std::string> fFilterLabels;
73 
78 
79  // future alg const std::string kKNNWithENu = "knnwithenu";
80  // future alg const std::string kKNNNoENu = "knnnoenu";
81 
83  int kIDNuE = 12; // TEMP, should find where stadard defs are.
84  /// Particle ID alorithm (loglikelihoods and dE/dx)
87 
88  /// MVAAlgorithm (Neural Net, KNN, etc. May include any algorithm that is
89  /// a subclass of slid::SPIDAlg)
92  /// FCL parameters for NuE Energy alorithm
94 
95  /// FCL parameters for particle ID alorithm (loglikelihoods and dE/dx)
97 
98  /// FCL parameters for SPID alorithm
100 
101  /// determine if library load is necessary or not when moving to a new run
102  bool libLoad = false;
103 
104  };
105  bool CompareByE( const slid::ShowerPID& a, const slid::ShowerPID& b);
106 } // namespace
107 
108 
109 
110 
111 namespace slid {
112  int Ncount0 = 0, Ncount1 = 0, Ncount2 = 0;
113  int Nfill = 0, Nput = 0;
114  //......................................................................
115 
117  : EDProducer(pset)
118  , fParticleIDAlg(0), fNuEEnergyAlg(0),
119  fMVAAlgEPi0(0),
120  fMVAAlgEPi0EL(0),
121  fNuEEnergyAlgPSet(pset.get< fhicl::ParameterSet >("NuEEnergyAlgPSet")),
122  fParticleIDAlgPSet(pset.get< fhicl::ParameterSet >("ParticleIDAlgPSet")),
123  fSPIDAlgPSet(pset.get< fhicl::ParameterSet >("SPIDAlgPSet"))
124  {
125  reconfigure(pset);
126  // Define output structures produced.
127  produces<std::vector<slid::ShowerPID> >();
128  produces< art::Assns<slid::ShowerPID, rb::Shower> >();
129  produces< art::Assns<slid::ShowerPID, rb::Cluster> >();
130  }
131 
132  //......................................................................
133 
135  {
136  if (fParticleIDAlg) delete fParticleIDAlg;
137  if (fNuEEnergyAlg) delete fNuEEnergyAlg;
138  if (fMVAAlgEPi0) delete fMVAAlgEPi0;
139  if (fMVAAlgEPi0EL) delete fMVAAlgEPi0EL;
140  }
141 
142  //......................................................................
143 
145  {
146  //
147  // Read FCL parameters
148  //
149  fNuEEnergyAlgPSet = pset.get< fhicl::ParameterSet >("NuEEnergyAlgPSet");
150  fParticleIDAlgPSet = pset.get< fhicl::ParameterSet >("ParticleIDAlgPSet");
151  fSPIDAlgPSet = pset.get< fhicl::ParameterSet >("SPIDAlgPSet");
152  fLibPath = util::EnvExpansion(pset.get< std::string >("LibraryPath"));
153  fSlicerLabel = pset.get< std::string >("SlicerLabel");
154  fShowerLabel = pset.get< std::string >("ShowerLabel");
155  fVertexLabel = pset.get< std::string >("VertexLabel");
156  fPIDOnlyForHighestEShower = pset.get<bool>("PIDOnlyForHighestEShower");
157  fObeyPreselection = pset.get<bool>("ObeyPreselection");
158  fSkipNoiseSlices = pset.get<bool>("SkipNoiseSlices");
159  fTrainingMode = pset.get<bool>("TrainingMode");
160  fFilterLabels = pset.get< std::vector<std::string> >("FilterLabels");
161  }
162 
164  {
166  novadaq::cnv::DetId detId = fGeom->DetId();
167 
168  //NOTE: This code is here because the geometry service only becomes available to determine the
169  //detector id at the begin run step. However at the present moment the libraries being loaded
170  //by the particle id and ann algorithms will not very from run to run and only depend on the detector id.
171  //If there are more then one run in a file this loading process takes a lot of time, only do it once per job.
172  if (!libLoad){
173  if (fParticleIDAlg) delete fParticleIDAlg;
174  if (fNuEEnergyAlg) delete fNuEEnergyAlg;
175  if (fMVAAlgEPi0) delete fMVAAlgEPi0;
176  if (fMVAAlgEPi0EL) delete fMVAAlgEPi0EL;
179 
180  if (!fTrainingMode){
183  }
184  libLoad = true;
185  }
186  }
187 
188 
189  //......................................................................
190 
192  {
193 
195 
196  MF_LOG_DEBUG("SPIDBuilderDEBUG") << "In SPIDBuilder::produce" << '\n';
197 
198  std::unique_ptr< std::vector<slid::ShowerPID> > showerpidcol(new std::vector<slid::ShowerPID>);
199  std::unique_ptr< art::Assns<slid::ShowerPID, rb::Shower> > assnSpidShower(new art::Assns<slid::ShowerPID, rb::Shower>);
200  std::unique_ptr< art::Assns<slid::ShowerPID, rb::Cluster> > assnSpidSlice(new art::Assns<slid::ShowerPID, rb::Cluster>);
201 
203  evt.getByLabel(fSlicerLabel, slicecol);
204  art::PtrVector<rb::Cluster> slicelist;
205  for (unsigned int i = 0; i < slicecol->size(); ++i) {
206  art::Ptr<rb::Cluster> slice(slicecol, i);
207  slicelist.push_back(slice);
208  }
209  //
210  // Grab list of rb::Shower* for calculations and reading and list of art::Ptr<rb::Shower>
211  // for writing the output data products. Since these **should** be indexed the same, this
212  // should be OK.
213  //
214  art::FindManyP<rb::Shower> fmtP(slicecol, evt, fShowerLabel);
216  art::FindManyP<rb::Vertex> fmv(slicecol, evt, fVertexLabel);
217 
218  // loop over all of the slices
219  for (unsigned int iSlice = 0; iSlice < slicelist.size(); ++iSlice) {
220  const size_t nextShwIdx = showerpidcol->size();
221  //
222  // Skip noise slices if this option is selected in FCL
223  //
224  if (fSkipNoiseSlices && slicelist[iSlice]->IsNoise()) continue;
225  //
226  // If there was nue preselection, skip slice if it failed if this option is selected in FCL
227  //
228  if (fObeyPreselection &&
229  rb::IsFiltered(evt, slicecol, iSlice, fFilterLabels))
230  continue;
231 
232  // get all of the showers in this slice
233  if ((!fmt.isValid())||(!fmtP.isValid())) continue;
234  const std::vector< const rb::Shower* > showerCollection = fmt.at(iSlice);
235  const std::vector < art::Ptr < rb::Shower >> showerCollectionPtr = fmtP.at(iSlice);
236  if ((showerCollection.empty())||(showerCollectionPtr.empty())) continue;
237 
238  //hold spid objects temporarily before sorting
239  std::vector<slid::ShowerPID> tempspidcol;
240  //keep track of shower index associationed with the spid object after sorting
241  std::vector<std::pair<double,int> > associationlist;
242 
243 
244  //get the vertex for the slice
245  const std::vector<art::Ptr<rb::Vertex>> vertCol = fmv.at(iSlice);
246  //currently only one vertex per slice
247  TVector3 vert = vertCol[0]->GetXYZ();
248 
249  //
250  // Get indices of highest and second-highest energy showers
251  //
252  double eShHighest = 0;
253  double sliceShwE = 0;
254  unsigned int iHighestEnergySh = -1;
255  for (unsigned int jShower = 0; jShower < showerCollection.size(); ++jShower) {
256  double thisESh = fNuEEnergyAlg->ShowerEnergy(showerCollection[jShower], showerCollection, evt);
257  sliceShwE += thisESh;
258  if (thisESh > eShHighest) {
259  eShHighest = thisESh;
260  iHighestEnergySh = jShower;
261  }
262  }
263 
264  TLorentzVector evtP4(0.0,0.0,0.0,0.0);
265  //
266  // Loop over showers to fill PID info
267  //
268  std::map<int, int> showerLidToShower;
269  unsigned int nShowers = showerCollection.size();
270  for (unsigned int iShower = 0; iShower < nShowers; ++iShower) {
271 
272 
273  slid::ShowerPID showerpid;
274 
275  double shwE = fNuEEnergyAlg->ShowerEnergy(showerCollection.at(iShower), showerCollection, evt);
276  if (std::isnan(sliceShwE) || std::isnan(shwE) || std::isinf(sliceShwE) || std::isinf(shwE)) {
277  mf::LogError("LIDBuilder")<< "sliceShwE/shwE is nan or inf";
278  continue;
279  }
280  if (shwE < 1e-10) shwE = 0.0;
281  showerpid.SetShowerEnergy(shwE);
282 
283  showerpid.SetDedx0(0);
284  showerpid.SetDedx1(0);
285  showerpid.SetDedx2(0);
286  showerpid.SetDedx3(0);
287  showerpid.SetDedx4(0);
288  showerpid.SetDedx5(0);
289  fParticleIDAlg->SetShower(showerCollection[iShower], showerCollection, evt);
290 
291  showerpid.SetELLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kELECTRON, showerCollection[iShower], showerCollection, evt));
292  showerpid.SetELLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kELECTRON, showerCollection[iShower], showerCollection, evt));
293  showerpid.SetGLLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kPHOTON, showerCollection[iShower], showerCollection, evt));
294  showerpid.SetGLLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kPHOTON, showerCollection[iShower], showerCollection, evt));
295  showerpid.SetMuLLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kMUON, showerCollection[iShower], showerCollection, evt));
296  showerpid.SetMuLLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kMUON, showerCollection[iShower],showerCollection, evt));
297  showerpid.SetPi0LLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kPI0, showerCollection[iShower],showerCollection, evt));
298  showerpid.SetPi0LLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kPI0, showerCollection[iShower],showerCollection, evt));
299  showerpid.SetPLLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kPROTON, showerCollection[iShower],showerCollection, evt));
300  showerpid.SetPLLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kPROTON, showerCollection[iShower], showerCollection, evt));
301  showerpid.SetNLLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kNEUTRON, showerCollection[iShower], showerCollection, evt));
302  showerpid.SetNLLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kNEUTRON, showerCollection[iShower],showerCollection, evt));
303  showerpid.SetPiLLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kPION, showerCollection[iShower],showerCollection, evt));
304  showerpid.SetPiLLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kPION, showerCollection[iShower],showerCollection, evt));
305 
306  if(showerCollection[iShower]->ExtentPlane()>0)showerpid.SetDedx0(fParticleIDAlg->PlaneLongDedx(0));
307  if(showerCollection[iShower]->ExtentPlane()>1)showerpid.SetDedx1(fParticleIDAlg->PlaneLongDedx(1));
308  if(showerCollection[iShower]->ExtentPlane()>2)showerpid.SetDedx2(fParticleIDAlg->PlaneLongDedx(2));
309  if(showerCollection[iShower]->ExtentPlane()>3)showerpid.SetDedx3(fParticleIDAlg->PlaneLongDedx(3));
310  if(showerCollection[iShower]->ExtentPlane()>4)showerpid.SetDedx4(fParticleIDAlg->PlaneLongDedx(4));
311  if(showerCollection[iShower]->ExtentPlane()>5)showerpid.SetDedx5(fParticleIDAlg->PlaneLongDedx(5));
312 //std::cout<<"SPID fParticleIDAlg dedx0 "<<fParticleIDAlg->PlaneLongDedx(0)<<std::endl;
313 //std::cout<<"SPID fParticleIDAlg ELLL, ELLT "<<showerpid.ELLL()<<" "<<showerpid.ELLT()<<std::endl;
314 //std::cout<<"SPID fParticleIDAlg GLLL, GLLT "<<showerpid.GLLL()<<" "<<showerpid.GLLT()<<std::endl;
315 
316 
317  TLorentzVector showerptrk(showerCollection[iShower]->Dir()[0]*shwE,
318  showerCollection[iShower]->Dir()[1]*shwE,
319  showerCollection[iShower]->Dir()[2]*shwE,
320  shwE);
321  TVector3 nuDir = fGeom->NuMIBeamDirection();
322  double theta = showerptrk.Angle(nuDir);
323  showerpid.SetCosTheta(cos(theta));
324 
325 
326  showerpid.SetVtxDoca(fParticleIDAlg->PointDoca(vertCol.at(0)->GetXYZ(),showerCollection[iShower]->Start(),showerCollection[iShower]->Stop()));
327 
328 
329  if (!fPIDOnlyForHighestEShower || iShower == iHighestEnergySh){
330  //in training mode just write -5, since the ann weights have not been determined
331  if (fTrainingMode){
332  showerpid.SetAlgDescription(""); // not needed?
333  showerpid.SetValEPi0(-5);
334  showerpid.SetValEPi0EL(-5);
335  }
336  //standard running, calculate the ann
337  else{
338  showerpid.SetValEPi0(fMVAAlgEPi0->CalcMVAResult(showerpid));
339  showerpid.SetValEPi0EL(fMVAAlgEPi0EL->CalcMVAResult(showerpid));
340  }
341  }
342 
343  showerpid.SetVtxDoca(fParticleIDAlg->PointDoca(vertCol.at(0)->GetXYZ(),showerCollection[iShower]->Start(),showerCollection[iShower]->Stop()));
344 
347  tempspidcol.push_back(showerpid);
348  associationlist.push_back(std::make_pair(showerpid.ShowerEnergy(),iShower));
349  }// end loop over showers
350 
351  std::sort( tempspidcol.begin(), tempspidcol.end(), CompareByE);
352  std::sort(associationlist.rbegin(), associationlist.rend());
353 
354  for (unsigned int ispid=0; ispid<tempspidcol.size(); ++ispid){
355  showerpidcol->push_back(tempspidcol[ispid]);
356  util::CreateAssn(evt, *(showerpidcol.get()), showerCollectionPtr[associationlist[ispid].second], *(assnSpidShower.get()));
357  }
358 
359  //if we have added new showers, make event pid and associations to slice
360  if (showerpidcol->size() > nextShwIdx){
361  //create new associations for the most recent set of showers
362  util::CreateAssn(evt, *showerpidcol, slicelist[iSlice], *(assnSpidSlice.get()), nextShwIdx, showerpidcol->size());
363  }
364 
365  } //end loop over slices
366 
367  evt.put(std::move(showerpidcol));
368  evt.put(std::move(assnSpidShower));
369  evt.put(std::move(assnSpidSlice));
370 
371  }// end of producer
372 
373  ////////////////////////////////////////////////////////////////////////
374 
376 
377  {
378  return a.ShowerEnergy() > b.ShowerEnergy();
379  }
380 
382 }
void SetPLLL(float in)
Definition: ShowerPID.h:48
void SetDedx5(float in)
Definition: ShowerPID.h:38
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.
fhicl::ParameterSet fSPIDAlgPSet
FCL parameters for SPID alorithm.
void SetVtxDoca(float in)
Definition: ShowerLID.h:84
bool libLoad
determine if library load is necessary or not when moving to a new run
This is a helper class for ParticleIDAlg that provides a tidy structure in which to hold the dE/dx hi...
pdg code and pid value
void SetDedx4(float in)
Definition: ShowerPID.h:37
double DedxTransLL(const slid::DedxParticleType partHyp, const rb::Shower *vShower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Calculate transverse dE/dx log-likelihood for specific particle hypothesis.
virtual void produce(art::Event &evt)
void SetDedx1(float in)
Definition: ShowerPID.h:34
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
std::map< int, float > fPartLongLL
Map of the longitudinal ll by paricle type.
Definition: ShowerLID.h:238
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
void SetPiLLT(float in)
Definition: ShowerPID.h:53
fhicl::ParameterSet fParticleIDAlgPSet
FCL parameters for particle ID alorithm (loglikelihoods and dE/dx)
std::vector< std::string > fFilterLabels
void SetShowerEnergy(float in)
Definition: ShowerLID.h:94
void SetDedx3(float in)
Definition: ShowerPID.h:36
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
bool CompareByE(const slid::ShowerLID &a, const slid::ShowerLID &b)
void SetValEPi0(float in)
Setters.
Definition: ShowerPID.h:30
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
void SetPi0LLL(float in)
Definition: ShowerPID.h:46
DEFINE_ART_MODULE(TestTMapFile)
double PlaneLongDedx(unsigned int pIdx)
return longitudinal dedx for a specified plane index
void SetNLLT(float in)
Definition: ShowerPID.h:51
slid::SPIDAlg * fMVAAlgEPi0
Definition: Run.h:21
void SetAlgDescription(const std::string &algDescription)
set algorithm description
Definition: ShowerLID.h:90
double CalcMVAResult(slid::ShowerPID &showerspid)
Calculate value of predictors for MVA Algorithm using candidate rb::Shower. A pointer to the object o...
Definition: SPIDAlg.cxx:247
ParticleIDAlg * fParticleIDAlg
Particle ID alorithm (loglikelihoods and dE/dx)
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
slid::SPIDAlg * fMVAAlgEPi0EL
void SetMuLLT(float in)
Definition: ShowerPID.h:45
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
const double a
T get(std::string const &key) const
Definition: ParameterSet.h:231
void SetPi0LLT(float in)
Definition: ShowerPID.h:47
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
std::map< int, float > fPartTransLL
Map of the transverse ll by particle type.
Definition: ShowerLID.h:241
void SetCosTheta(float in)
Definition: ShowerLID.h:47
fhicl::ParameterSet fNuEEnergyAlgPSet
FCL parameters for NuE Energy alorithm.
std::map< int, float > fPartTransLL
Map of the transverse ll by particle type.
Vertex location in position and time.
double PointDoca(TVector3 vtx, TVector3 start, TVector3 stop)
Get distance of closest approach between shower and vertex.
SPIDBuilder(const fhicl::ParameterSet &pset)
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:302
void SetPiLLL(float in)
Definition: ShowerPID.h:52
double ShowerEnergy(const rb::Shower *shower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Returns the total energy of a shower. along with corrections due to dead material and threshold effec...
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
void SetDedx2(float in)
Definition: ShowerPID.h:35
NuEEnergyAlg * fNuEEnergyAlg
Calculate dE/dx and log likelihoods for different particle hypotheses. This information will be of us...
float ShowerEnergy() const
Definition: ShowerLID.h:157
Calculates deposited and corrected energy of the electron shower and of electron flavoured neutrino...
void SetMuLLL(float in)
Definition: ShowerPID.h:44
#define MF_LOG_DEBUG(id)
void SetELLL(float in)
Definition: ShowerPID.h:40
const hit & b
Definition: hits.cxx:21
void SetValEPi0EL(float in)
Definition: ShowerPID.h:31
void SetGLLL(float in)
Definition: ShowerPID.h:42
T cos(T number)
Definition: d0nt_math.hpp:78
Build slid::LID objects to store electron ID, if asked for, otherwise, calculate LID info and make av...
Definition: FillPIDs.h:13
virtual void beginRun(art::Run &run)
Float_t e
Definition: plot.C:35
std::map< int, float > fPartLongLL
Map of the longitudinal ll by paricle type.
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
void SetPLLT(float in)
Definition: ShowerPID.h:49
void SetELLT(float in)
Definition: ShowerPID.h:41
void SetGLLT(float in)
Definition: ShowerPID.h:43
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual void reconfigure(const fhicl::ParameterSet &pset)
void SetNLLL(float in)
Definition: ShowerPID.h:50
void SetDedx0(float in)
Definition: ShowerPID.h:33
double DedxLongLL(const slid::DedxParticleType partHyp, const rb::Shower *vShower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis.
void SetShower(const rb::Shower *vShower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Set rb::Prong to be analyzed. This must be set before any calculations can be done.
enum BeamMode string