MakeNCPi0BkgRej_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file MakeNCPi0BkgRej_module.cc
3 // \brief A module to write out the NCPi0 object
4 //
5 // \author Daisy Kalra - kalra@fnal.gov
6 ////////////////////////////////////////////////////////////////////////
7 // art includes
18 #include "cetlib/search_path.h"
19 #include "cetlib/filesystem.h"
20 #include "fhiclcpp/ParameterSet.h"
22 #include "cetlib_except/exception.h"
23 #include <memory>
24 #include "TMath.h"
25 #include "TVector3.h"
26 // TMVA-Based
27 #include "TMVA/Reader.h"
28 
29 // novasoft includes
31 #include "CVN/func/Result.h"
33 #include "ShowerLID/ShowerLID.h"
34 #include "Geometry/Geometry.h"
36 #include "Utilities/AssociationUtil.h"
37 #include "RecoBase/Cluster.h"
38 #include "RecoBase/Shower.h"
39 #include "RecoBase/Vertex.h"
40 #include "RecoBase/FilterList.h"
42 #include "NueSandbox/NueSandObj.h"
43 
44 
45 namespace ncpi0
46 {
48  public:
49 
50  explicit MakeNCPi0BkgRej(fhicl::ParameterSet const &p);
52 
53  void beginRun(art::Run &run) override;
54  void reconfigure(const fhicl::ParameterSet &p);
55  void produce(art::Event &evt) override;
56 
57  private:
58 
69  bool fObeyPreselection; // nue presel: no
70  std::vector<std::string> fPreselectionLabels; //veto: yes
71 
72 
73  TMVA::Reader* fReaderBDTG;
74  TMVA::Reader* fReaderBDTGLowThres;
75 
76  // Variables need to send to TMVA::Reader
77  float TMVAvarsBDTG[9];
79  };
80 }
81 
82 //////////////////////////////////////////////////////////////////////////
83 namespace ncpi0
84 {
85 
86  //////////////////////////////////////////////////////////////////////////
88  : fReaderBDTG(nullptr) //for NCPi0
89  {
90  this->reconfigure(p);
91  this->produces< std::vector<ncpi0::NCPi0BkgRej> > ();
92  this->produces< art::Assns<ncpi0::NCPi0BkgRej, rb::Cluster> > ();
93  }
94 
95  //////////////////////////////////////////////////////////////////////////
97  {
98  if(fReaderBDTG) delete fReaderBDTG;
99  }
100 
101  //////////////////////////////////////////////////////////////////////////
103  {
104 
105  fReaderBDTG = new TMVA::Reader;
106  fReaderBDTG->AddVariable("cvnnumu", &TMVAvarsBDTG[0]);
107  fReaderBDTG->AddVariable("prong1epi0LLL", &TMVAvarsBDTG[1]);
108  fReaderBDTG->AddVariable("prong1epiLLL", &TMVAvarsBDTG[2]);
109  fReaderBDTG->AddVariable("prong1ContPl", &TMVAvarsBDTG[3]);
110  fReaderBDTG->AddVariable("prong1epLLT", &TMVAvarsBDTG[4]);
111  fReaderBDTG->AddVariable("prong1Width", &TMVAvarsBDTG[5]);
112  fReaderBDTG->AddVariable("cvnncid", &TMVAvarsBDTG[6]);
113  fReaderBDTG->AddVariable("prong2dedx", &TMVAvarsBDTG[7]);
114  fReaderBDTG->AddVariable("prong1MissingPl", &TMVAvarsBDTG[8]);
115 
116 
118  // Make sure we can find the Gradient Boosted Decision (For NCPi0) Trees
119  // weight file before we set up MVA
120  if(!cet::file_exists(pidlibpath+fNCPi0BkgRejPIDFile)){
121  throw cet::exception("MakeNCPi0BkgRej")
122  << "Couldn't find file " << pidlibpath+fNCPi0BkgRejPIDFile <<std::endl;
123  }
124  fReaderBDTG->BookMVA("BDTG", pidlibpath+fNCPi0BkgRejPIDFile);
125 
126  //For 10 variables PID with 100meV latest weight file
127 
128  fReaderBDTGLowThres = new TMVA::Reader;
129  fReaderBDTGLowThres->AddVariable("cvnnumu", &TMVAvarsBDTG[0]);
130  fReaderBDTGLowThres->AddVariable("prong1epi0LLL", &TMVAvarsBDTG[1]);
131  fReaderBDTGLowThres->AddVariable("prong1epiLLL", &TMVAvarsBDTG[2]);
132  fReaderBDTGLowThres->AddVariable("prong1ContPl", &TMVAvarsBDTG[3]);
133  fReaderBDTGLowThres->AddVariable("prong1epLLT", &TMVAvarsBDTG[4]);
134  fReaderBDTGLowThres->AddVariable("prong1Width", &TMVAvarsBDTG[5]);
135  fReaderBDTGLowThres->AddVariable("cvnncid", &TMVAvarsBDTG[6]);
136  fReaderBDTGLowThres->AddVariable("prong1dedx", &TMVAvarsBDTG[7]);
137  fReaderBDTGLowThres->AddVariable("prong2dedx", &TMVAvarsBDTG[8]);
138  fReaderBDTGLowThres->AddVariable("prong1MissingPl", &TMVAvarsBDTG[9]);
139 
140 
142  // Make sure we can find the Gradient Boosted Decision (For NCPi0) Trees
143  // weight file before we set up MVA
144  if(!cet::file_exists(pidlibpathLT+fNCPi0BkgRejPIDFileLT)){
145  throw cet::exception("MakeNCPi0BkgRejLT")
146  << "Couldn't find file " << pidlibpathLT+fNCPi0BkgRejPIDFileLT <<std::endl;
147  }
148  fReaderBDTGLowThres->BookMVA("BDTG", pidlibpathLT+fNCPi0BkgRejPIDFileLT);
149 
150 
151  }
152 
153  //////////////////////////////////////////////////////////////////////////
155  {
156  fSliceLabel = p.get< std::string >("SliceLabel");
157  fCVNLabel = p.get< std::string >("CVNLabel");
158  fShowerLabel = p.get< std::string >("ShowerLabel");
159  fShowerLIDLabel = p.get< std::string >("ShowerLIDLabel");
160  fNueSandboxLabel = p.get< std::string >("NueSandboxLabel");
161  fPIDLibPath = p.get< std::string >("PIDLibPath");
162  fNCPi0BkgRejPIDFile = p.get< std::string >("NCPi0BkgRejPIDFile");
163  fPIDLibPathLT = p.get< std::string >("PIDLibPathLT");
164  fNCPi0BkgRejPIDFileLT = p.get< std::string >("NCPi0BkgRejPIDFileLT");
165  fObeyPreselection = p.get< bool >("ObeyPreselection");
166  fPreselectionLabels = p.get< std::vector<std::string> > ("PreselectionLabels");
167  }
168 
169  //////////////////////////////////////////////////////////////////////////
171  {
172  float prong1dedx = -5.;
173  float prong2dedx = -5.;
174  float prong1epi0LLL = -5.;
175  float prong1epLLT = -5.;
176  float prong1epiLLL = -5.;
177  float prong1ContPl = -5.;
178  float prong1MissingPl = -5.;
179  float prong1Width = -5.;
180 
181  // Get the slices in the event
183  evt.getByLabel( fSliceLabel, sliceHandle);
184  std::vector<art::Ptr<rb::Cluster> > slices;
185  art::fill_ptr_vector(slices, sliceHandle);
186 
187  if(sliceHandle.failedToGet()){
188  mf::LogInfo("MakeNCPi0BkgRej") << "NO SLICE FOUND, skipping event.....";
189  }
190 
191 
192  std::unique_ptr< std::vector<ncpi0::NCPi0BkgRej> > ncpi0BkgRejjies(new std::vector<ncpi0::NCPi0BkgRej>);
193  std::unique_ptr< art::Assns<ncpi0::NCPi0BkgRej, rb::Cluster> > sliceAssn(new art::Assns<ncpi0::NCPi0BkgRej, rb::Cluster>);
194 
195  ncpi0::NCPi0BkgRej ncpi0bkgrej;
196 
197  // Get things associated with slices
198  art::FindMany < slid::ShowerLID > showerLidAssn(sliceHandle, evt, fShowerLIDLabel);
199  unsigned int nslices = slices.size();
200 
201  for(unsigned int isli = 0; isli < nslices; ++isli){
202  if(fObeyPreselection && rb::IsFiltered(evt, sliceHandle, isli, fPreselectionLabels)) continue;
203  //if(fObeyPreselection && rb::IsFiltered(evt, sliceHandle, isli)) continue;
204  if(slices.empty()) continue;
205 
206  art::Ptr<rb::Cluster> thisSlice = slices.at(isli);
207  if(thisSlice->IsNoise()) continue;
208 
209  // Only execute when there are showers associated with the shwlids
210  if(!showerLidAssn.isValid()) continue;
211  std::vector< const slid::ShowerLID* > shwlids = showerLidAssn.at(isli);
212  if(shwlids.empty()) continue;
213 
214  sort(shwlids.begin(), shwlids.end(), slid::CompareByEnergy);
215  art::FindOneP<rb::Shower> foSh(shwlids, evt, fShowerLIDLabel);
216  if(!foSh.isValid()) continue;
217 
218  //Getting 3dProng1 variables
219  std::vector< art::Ptr<rb::Shower> > showers;
220  cet::maybe_ref< art::Ptr<rb::Shower> const > roShws(foSh.at(0));
221  art::Ptr<rb::Shower> showers_dist = roShws.ref();
222  showers.push_back(showers_dist);
223  art::FindOneP<slid::ShowerLID> foShLID(showers, evt,fShowerLIDLabel);
224  if(!foShLID.isValid()) continue;
225  cet::maybe_ref< art::Ptr<slid::ShowerLID> const > roShLID(foShLID.at(0));
226  art::Ptr<slid::ShowerLID> showerLID = roShLID.ref();
227 
228  prong1epi0LLL = showerLID->EPi0LLL();
229  prong1epLLT = showerLID->EPLLT();
230  prong1epiLLL = showerLID->EPiLLL();
231  prong1ContPl = showers_dist->MostContiguousPlanes(geo::kXorY);
232  prong1MissingPl = showers_dist->MostMissingPlanes(geo::kXorY);
233  prong1Width = showerLID->Radius();
234 
235  //Prong1 and Prong2 dEdx as a input to Grad BDT
236 
237  //
238  // ATTENTION:
239  //
240  // Warning future users about using the nue sandbox...
241  std::cout << "\n\n\nWARNING: This module is attempting to use the nue sandbox variables which are no longer a part of standard production.\n\n\n";
242 
243  art::FindOneP<nuesand::NueSandObj> fmNueSand(sliceHandle, evt, fNueSandboxLabel);
244  if(!fmNueSand.isValid()) continue;
245  art::Ptr<nuesand::NueSandObj> nuesands = fmNueSand.at(isli);
246  if(!nuesands) continue;
247 
248  prong1dedx = nuesands->fdEdxProng1;
249  prong2dedx = nuesands->fdEdxProng2;
250 
251  // Need CVN, input to Grad BDT
252 
253  art::FindOneP<cvn::Result> fmCVN(sliceHandle, evt, fCVNLabel);
254  if(!fmCVN.isValid()) continue; // CVN Validation
255  art::Ptr<cvn::Result> cvns = fmCVN.at(isli);
256  if(!cvns) continue;
257  float cvnnumu = 0.;
258  float cvnncid = 0.;
259 
260  //
261  // ATTENTION:
262  //
263  // Warning future users about using CVN with the wrong number of output labels...
264  std::cout << "\n\n\nWARNING: This module is assuming an OLD output label scheme for CVN. The number of output labels has not been 14 since prod3. Any results obtained this way are likely nonsense.\n\n\n";
265 
266  // CVN values
267  cvnnumu = cvns->fOutput[0] + cvns->fOutput[1] + cvns->fOutput[2] + cvns->fOutput[3];
268  cvnncid = cvns->fOutput[13];
269 
270  //PID with K.E threshold > 0.5 GeV and with Reco K.E Cut > 0.5GeV
271  TMVAvarsBDTG[0] = cvnnumu;
272  TMVAvarsBDTG[1] = prong1epi0LLL;
273  TMVAvarsBDTG[2] = prong1epiLLL;
274  TMVAvarsBDTG[3] = prong1ContPl;
275  TMVAvarsBDTG[4] = prong1epLLT;
276  TMVAvarsBDTG[5] = prong1Width;
277  TMVAvarsBDTG[6] = cvnncid;
278  TMVAvarsBDTG[7] = prong2dedx;
279  TMVAvarsBDTG[8] = prong1MissingPl;
280  ncpi0bkgrej.SetNCPIDBDTG(fReaderBDTG->EvaluateMVA("BDTG"));
281 
282  //PID with lower threshold > 0.2 GeV and no reco K.E cut
283  TMVAvarsBDTGLowThres[0] = cvnnumu;
284  TMVAvarsBDTGLowThres[1] = prong1epi0LLL;
285  TMVAvarsBDTGLowThres[2] = prong1epiLLL;
286  TMVAvarsBDTGLowThres[3] = prong1ContPl;
287  TMVAvarsBDTGLowThres[4] = prong1epLLT;
288  TMVAvarsBDTGLowThres[5] = prong1Width;
289  TMVAvarsBDTGLowThres[6] = cvnncid;
290  TMVAvarsBDTGLowThres[7] = prong1dedx;
291  TMVAvarsBDTGLowThres[8] = prong2dedx;
292  TMVAvarsBDTGLowThres[9] = prong1MissingPl;
293  ncpi0bkgrej.SetNCPIDBDTGLT(fReaderBDTGLowThres->EvaluateMVA("BDTG"));
294 
295  ncpi0BkgRejjies->push_back(ncpi0bkgrej);
296 
297  util::CreateAssn(*this, evt, *ncpi0BkgRejjies, thisSlice ,*sliceAssn);
298 
299  } // end loop over slices
300 
301  evt.put(std::move(ncpi0BkgRejjies));
302  evt.put(std::move(sliceAssn));
303 
304  } // end of producer
305 
306 
307  ////////////////////////////////////////////////////////////////////////
308 
309 } // end namespace ncpi0
310 
312 
313 ////////////////////////////////////////////////////////////////////////
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.
X or Y views.
Definition: PlaneGeo.h:30
void SetNCPIDBDTGLT(double input)
Definition: NCPi0BkgRej.h:27
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
const char * p
Definition: xmltok.h:285
std::vector< std::string > fPreselectionLabels
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
DEFINE_ART_MODULE(TestTMapFile)
std::vector< float > fOutput
Vector of outputs from neural net.
Definition: Result.h:30
bool CompareByEnergy(const slid::ShowerLID *a, const slid::ShowerLID *b)
Definition: ShowerLID.cxx:51
void beginRun(art::Run &run) override
Definition: Run.h:31
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
Result for CVN.
T get(std::string const &key) const
Definition: ParameterSet.h:231
void reconfigure(const fhicl::ParameterSet &p)
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
void produce(art::Event &evt) override
Vertex location in position and time.
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
void SetNCPIDBDTG(double input)
Definition: NCPi0BkgRej.h:24
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Calculate dE/dx and log likelihoods for different particle hypotheses. This information will be of us...
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
MakeNCPi0BkgRej(fhicl::ParameterSet const &p)
bool file_exists(std::string const &qualified_filename)
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
Definition: FillPIDs.h:19
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196