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
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  : fParticleIDAlg(0), fNuEEnergyAlg(0),
118  fMVAAlgEPi0(0),
119  fMVAAlgEPi0EL(0),
120  fNuEEnergyAlgPSet(pset.get< fhicl::ParameterSet >("NuEEnergyAlgPSet")),
121  fParticleIDAlgPSet(pset.get< fhicl::ParameterSet >("ParticleIDAlgPSet")),
122  fSPIDAlgPSet(pset.get< fhicl::ParameterSet >("SPIDAlgPSet"))
123  {
124  reconfigure(pset);
125  // Define output structures produced.
126  produces<std::vector<slid::ShowerPID> >();
127  produces< art::Assns<slid::ShowerPID, rb::Shower> >();
128  produces< art::Assns<slid::ShowerPID, rb::Cluster> >();
129  }
130 
131  //......................................................................
132 
134  {
135  if (fParticleIDAlg) delete fParticleIDAlg;
136  if (fNuEEnergyAlg) delete fNuEEnergyAlg;
137  if (fMVAAlgEPi0) delete fMVAAlgEPi0;
138  if (fMVAAlgEPi0EL) delete fMVAAlgEPi0EL;
139  }
140 
141  //......................................................................
142 
144  {
145  //
146  // Read FCL parameters
147  //
148  fNuEEnergyAlgPSet = pset.get< fhicl::ParameterSet >("NuEEnergyAlgPSet");
149  fParticleIDAlgPSet = pset.get< fhicl::ParameterSet >("ParticleIDAlgPSet");
150  fSPIDAlgPSet = pset.get< fhicl::ParameterSet >("SPIDAlgPSet");
151  fLibPath = util::EnvExpansion(pset.get< std::string >("LibraryPath"));
152  fSlicerLabel = pset.get< std::string >("SlicerLabel");
153  fShowerLabel = pset.get< std::string >("ShowerLabel");
154  fVertexLabel = pset.get< std::string >("VertexLabel");
155  fPIDOnlyForHighestEShower = pset.get<bool>("PIDOnlyForHighestEShower");
156  fObeyPreselection = pset.get<bool>("ObeyPreselection");
157  fSkipNoiseSlices = pset.get<bool>("SkipNoiseSlices");
158  fTrainingMode = pset.get<bool>("TrainingMode");
159  fFilterLabels = pset.get< std::vector<std::string> >("FilterLabels");
160  }
161 
163  {
165  novadaq::cnv::DetId detId = fGeom->DetId();
166 
167  //NOTE: This code is here because the geometry service only becomes available to determine the
168  //detector id at the begin run step. However at the present moment the libraries being loaded
169  //by the particle id and ann algorithms will not very from run to run and only depend on the detector id.
170  //If there are more then one run in a file this loading process takes a lot of time, only do it once per job.
171  if (!libLoad){
172  if (fParticleIDAlg) delete fParticleIDAlg;
173  if (fNuEEnergyAlg) delete fNuEEnergyAlg;
174  if (fMVAAlgEPi0) delete fMVAAlgEPi0;
175  if (fMVAAlgEPi0EL) delete fMVAAlgEPi0EL;
178 
179  if (!fTrainingMode){
182  }
183  libLoad = true;
184  }
185  }
186 
187 
188  //......................................................................
189 
191  {
192 
194 
195  LOG_DEBUG("SPIDBuilderDEBUG") << "In SPIDBuilder::produce" << '\n';
196 
197  std::unique_ptr< std::vector<slid::ShowerPID> > showerpidcol(new std::vector<slid::ShowerPID>);
198  std::unique_ptr< art::Assns<slid::ShowerPID, rb::Shower> > assnSpidShower(new art::Assns<slid::ShowerPID, rb::Shower>);
199  std::unique_ptr< art::Assns<slid::ShowerPID, rb::Cluster> > assnSpidSlice(new art::Assns<slid::ShowerPID, rb::Cluster>);
200 
202  evt.getByLabel(fSlicerLabel, slicecol);
203  art::PtrVector<rb::Cluster> slicelist;
204  for (unsigned int i = 0; i < slicecol->size(); ++i) {
205  art::Ptr<rb::Cluster> slice(slicecol, i);
206  slicelist.push_back(slice);
207  }
208  //
209  // Grab list of rb::Shower* for calculations and reading and list of art::Ptr<rb::Shower>
210  // for writing the output data products. Since these **should** be indexed the same, this
211  // should be OK.
212  //
213  art::FindManyP<rb::Shower> fmtP(slicecol, evt, fShowerLabel);
215  art::FindManyP<rb::Vertex> fmv(slicecol, evt, fVertexLabel);
216 
217  // loop over all of the slices
218  for (unsigned int iSlice = 0; iSlice < slicelist.size(); ++iSlice) {
219  const size_t nextShwIdx = showerpidcol->size();
220  //
221  // Skip noise slices if this option is selected in FCL
222  //
223  if (fSkipNoiseSlices && slicelist[iSlice]->IsNoise()) continue;
224  //
225  // If there was nue preselection, skip slice if it failed if this option is selected in FCL
226  //
227  if (fObeyPreselection &&
228  rb::IsFiltered(evt, slicecol, iSlice, fFilterLabels))
229  continue;
230 
231  // get all of the showers in this slice
232  if ((!fmt.isValid())||(!fmtP.isValid())) continue;
233  const std::vector< const rb::Shower* > showerCollection = fmt.at(iSlice);
234  const std::vector < art::Ptr < rb::Shower >> showerCollectionPtr = fmtP.at(iSlice);
235  if ((showerCollection.empty())||(showerCollectionPtr.empty())) continue;
236 
237  //hold spid objects temporarily before sorting
238  std::vector<slid::ShowerPID> tempspidcol;
239  //keep track of shower index associationed with the spid object after sorting
240  std::vector<std::pair<double,int> > associationlist;
241 
242 
243  //get the vertex for the slice
244  const std::vector<art::Ptr<rb::Vertex>> vertCol = fmv.at(iSlice);
245  //currently only one vertex per slice
246  TVector3 vert = vertCol[0]->GetXYZ();
247 
248  //
249  // Get indices of highest and second-highest energy showers
250  //
251  double eShHighest = 0;
252  double sliceShwE = 0;
253  unsigned int iHighestEnergySh = -1;
254  for (unsigned int jShower = 0; jShower < showerCollection.size(); ++jShower) {
255  double thisESh = fNuEEnergyAlg->ShowerEnergy(showerCollection[jShower], showerCollection, evt);
256  sliceShwE += thisESh;
257  if (thisESh > eShHighest) {
258  eShHighest = thisESh;
259  iHighestEnergySh = jShower;
260  }
261  }
262 
263  TLorentzVector evtP4(0.0,0.0,0.0,0.0);
264  //
265  // Loop over showers to fill PID info
266  //
267  std::map<int, int> showerLidToShower;
268  unsigned int nShowers = showerCollection.size();
269  for (unsigned int iShower = 0; iShower < nShowers; ++iShower) {
270 
271 
272  slid::ShowerPID showerpid;
273 
274  double shwE = fNuEEnergyAlg->ShowerEnergy(showerCollection.at(iShower), showerCollection, evt);
275  if (std::isnan(sliceShwE) || std::isnan(shwE) || std::isinf(sliceShwE) || std::isinf(shwE)) {
276  mf::LogError("LIDBuilder")<< "sliceShwE/shwE is nan or inf";
277  continue;
278  }
279  if (shwE < 1e-10) shwE = 0.0;
280  showerpid.SetShowerEnergy(shwE);
281 
282  showerpid.SetDedx0(0);
283  showerpid.SetDedx1(0);
284  showerpid.SetDedx2(0);
285  showerpid.SetDedx3(0);
286  showerpid.SetDedx4(0);
287  showerpid.SetDedx5(0);
288  fParticleIDAlg->SetShower(showerCollection[iShower], showerCollection, evt);
289 
290  showerpid.SetELLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kELECTRON, showerCollection[iShower], showerCollection, evt));
291  showerpid.SetELLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kELECTRON, showerCollection[iShower], showerCollection, evt));
292  showerpid.SetGLLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kPHOTON, showerCollection[iShower], showerCollection, evt));
293  showerpid.SetGLLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kPHOTON, showerCollection[iShower], showerCollection, evt));
294  showerpid.SetMuLLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kMUON, showerCollection[iShower], showerCollection, evt));
295  showerpid.SetMuLLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kMUON, showerCollection[iShower],showerCollection, evt));
296  showerpid.SetPi0LLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kPI0, showerCollection[iShower],showerCollection, evt));
297  showerpid.SetPi0LLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kPI0, showerCollection[iShower],showerCollection, evt));
298  showerpid.SetPLLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kPROTON, showerCollection[iShower],showerCollection, evt));
299  showerpid.SetPLLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kPROTON, showerCollection[iShower], showerCollection, evt));
300  showerpid.SetNLLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kNEUTRON, showerCollection[iShower], showerCollection, evt));
301  showerpid.SetNLLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kNEUTRON, showerCollection[iShower],showerCollection, evt));
302  showerpid.SetPiLLL(fParticleIDAlg->DedxLongLL(slid::DedxParticleType::kPION, showerCollection[iShower],showerCollection, evt));
303  showerpid.SetPiLLT(fParticleIDAlg->DedxTransLL(slid::DedxParticleType::kPION, showerCollection[iShower],showerCollection, evt));
304 
305  if(showerCollection[iShower]->ExtentPlane()>0)showerpid.SetDedx0(fParticleIDAlg->PlaneLongDedx(0));
306  if(showerCollection[iShower]->ExtentPlane()>1)showerpid.SetDedx1(fParticleIDAlg->PlaneLongDedx(1));
307  if(showerCollection[iShower]->ExtentPlane()>2)showerpid.SetDedx2(fParticleIDAlg->PlaneLongDedx(2));
308  if(showerCollection[iShower]->ExtentPlane()>3)showerpid.SetDedx3(fParticleIDAlg->PlaneLongDedx(3));
309  if(showerCollection[iShower]->ExtentPlane()>4)showerpid.SetDedx4(fParticleIDAlg->PlaneLongDedx(4));
310  if(showerCollection[iShower]->ExtentPlane()>5)showerpid.SetDedx5(fParticleIDAlg->PlaneLongDedx(5));
311 //std::cout<<"SPID fParticleIDAlg dedx0 "<<fParticleIDAlg->PlaneLongDedx(0)<<std::endl;
312 //std::cout<<"SPID fParticleIDAlg ELLL, ELLT "<<showerpid.ELLL()<<" "<<showerpid.ELLT()<<std::endl;
313 //std::cout<<"SPID fParticleIDAlg GLLL, GLLT "<<showerpid.GLLL()<<" "<<showerpid.GLLT()<<std::endl;
314 
315 
316  TLorentzVector showerptrk(showerCollection[iShower]->Dir()[0]*shwE,
317  showerCollection[iShower]->Dir()[1]*shwE,
318  showerCollection[iShower]->Dir()[2]*shwE,
319  shwE);
320  TVector3 nuDir = fGeom->NuMIBeamDirection();
321  double theta = showerptrk.Angle(nuDir);
322  showerpid.SetCosTheta(cos(theta));
323 
324 
325  showerpid.SetVtxDoca(fParticleIDAlg->PointDoca(vertCol.at(0)->GetXYZ(),showerCollection[iShower]->Start(),showerCollection[iShower]->Stop()));
326 
327 
328  if (!fPIDOnlyForHighestEShower || iShower == iHighestEnergySh){
329  //in training mode just write -5, since the ann weights have not been determined
330  if (fTrainingMode){
331  showerpid.SetAlgDescription(""); // not needed?
332  showerpid.SetValEPi0(-5);
333  showerpid.SetValEPi0EL(-5);
334  }
335  //standard running, calculate the ann
336  else{
337  showerpid.SetValEPi0(fMVAAlgEPi0->CalcMVAResult(showerpid));
338  showerpid.SetValEPi0EL(fMVAAlgEPi0EL->CalcMVAResult(showerpid));
339  }
340  }
341 
342  showerpid.SetVtxDoca(fParticleIDAlg->PointDoca(vertCol.at(0)->GetXYZ(),showerCollection[iShower]->Start(),showerCollection[iShower]->Stop()));
343 
346  tempspidcol.push_back(showerpid);
347  associationlist.push_back(std::make_pair(showerpid.ShowerEnergy(),iShower));
348  }// end loop over showers
349 
350  std::sort( tempspidcol.begin(), tempspidcol.end(), CompareByE);
351  std::sort(associationlist.rbegin(), associationlist.rend());
352 
353  for (unsigned int ispid=0; ispid<tempspidcol.size(); ++ispid){
354  showerpidcol->push_back(tempspidcol[ispid]);
355  util::CreateAssn(*this, evt, *(showerpidcol.get()), showerCollectionPtr[associationlist[ispid].second], *(assnSpidShower.get()));
356  }
357 
358  //if we have added new showers, make event pid and associations to slice
359  if (showerpidcol->size() > nextShwIdx){
360  //create new associations for the most recent set of showers
361  util::CreateAssn(*this, evt, *showerpidcol, slicelist[iSlice], *(assnSpidSlice.get()), nextShwIdx, showerpidcol->size());
362  }
363 
364  } //end loop over slices
365 
366  evt.put(std::move(showerpidcol));
367  evt.put(std::move(assnSpidShower));
368  evt.put(std::move(assnSpidSlice));
369 
370  }// end of producer
371 
372  ////////////////////////////////////////////////////////////////////////
373 
375 
376  {
377  return a.ShowerEnergy() > b.ShowerEnergy();
378  }
379 
381 }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
void SetPLLL(float in)
Definition: ShowerPID.h:48
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
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
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:31
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
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
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
ParticleIDAlg * fParticleIDAlg
Particle ID alorithm (loglikelihoods and dE/dx)
slid::SPIDAlg * fMVAAlgEPi0EL
void SetMuLLT(float in)
Definition: ShowerPID.h:45
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
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:308
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
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
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
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Calculates deposited and corrected energy of the electron shower and of electron flavoured neutrino...
void SetMuLLL(float in)
Definition: ShowerPID.h:44
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.
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.