CVNEvaluator_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file CVNEvaluator_module.cc
3 // \brief Producer module creating CVN neural net results
4 // \author Alexander Radovic - a.radovic@gmail.com
5 // Saul Alonso Monsalve - saul.alonso.monsalve@cern.ch
6 ////////////////////////////////////////////////////////////////////////
7 
8 // C/C++ includes
9 #include <iostream>
10 #include <sstream>
11 
12 // Framework includes
16 #include "art_root_io/TFileDirectory.h"
17 #include "art_root_io/TFileService.h"
19 #include "fhiclcpp/ParameterSet.h"
23 
26 //#include "dunereco/CVN/art/CaffeNetHandler.h"
30 
31 namespace cvn {
32 
33  class CVNEvaluator : public art::EDProducer {
34  public:
35  explicit CVNEvaluator(fhicl::ParameterSet const& pset);
36  ~CVNEvaluator();
37 
38  void produce(art::Event& evt);
39  void beginJob();
40  void endJob();
41 
42 
43 
44  private:
45 
46  /// Module label for input pixel maps
49 
50  /// Can use Caffe or Tensorflow
52 
53  //cvn::CaffeNetHandler fCaffeHandler;
55 
56  /// Number of outputs fron neural net
57  //unsigned int fNOutput;
58 
59  /// If there are multiple pixel maps per event can we use them?
61 
62  unsigned int fTotal;
63  unsigned int fCorrect;
64  unsigned int fFullyCorrect;
65 
66  unsigned int fTotNue;
67  unsigned int fSelNue;
68  unsigned int fTotNumu;
69  unsigned int fSelNumu;
70 
71  // Three binned versions of above
72  std::vector<unsigned int> fTotNueBins;
73  std::vector<unsigned int> fSelNueBins;
74  std::vector<unsigned int> fTotNumuBins;
75  std::vector<unsigned int> fSelNumuBins;
76  };
77 
78  //.......................................................................
80  fPixelMapInput (pset.get<std::string> ("PixelMapInput")),
81  fResultLabel (pset.get<std::string> ("ResultLabel")),
82  fCVNType (pset.get<std::string> ("CVNType")),
83  //fCaffeHandler (pset.get<fhicl::ParameterSet> ("CaffeNetHandler")),
84  fTFHandler (pset.get<fhicl::ParameterSet> ("TFNetHandler")),
85  //fNOutput (fCaffeHandler.NOutput()),
86  fMultiplePMs (pset.get<bool> ("MultiplePMs"))
87  {
88  produces< std::vector<cvn::Result> >(fResultLabel);
89  fTotal = 0;
90  fCorrect = 0;
91  fFullyCorrect = 0;
92  fTotNue = 0;
93  fSelNue = 0;
94  fTotNumu = 0;
95  fSelNumu = 0;
96 
97  fTotNueBins = {0,0,0};
98  fSelNueBins = {0,0,0};
99 
100  fTotNumuBins = {0,0,0};
101  fSelNumuBins = {0,0,0};
102  }
103  //......................................................................
105  {
106  //======================================================================
107  // Clean up any memory allocated by your module
108  //======================================================================
109  }
110 
111  //......................................................................
113  { }
114 
115  //......................................................................
117  {
118  /*
119  float tot = static_cast<float>(fTotal);
120  std::cout << "Total of " << fTotal << " events passed the fiducial volume cut and were classified" << std::endl;
121  if(fTotal > 0) std::cout << "Correct flavour (interaction) identification = " << 100.*(fCorrect / tot) << "% (" << 100.*(fFullyCorrect/tot) << "%)" << std::endl;
122  if(fTotNue > 0) std::cout << "Nue efficiency = " << 100.*(fSelNue/static_cast<float>(fTotNue)) << "%" << std::endl;
123  for(unsigned int i = 0; i < fTotNueBins.size(); ++i){
124  if(fTotNueBins[i] > 0) std::cout << "Nue efficiency " << i << "= " << 100.*(fSelNueBins[i]/static_cast<float>(fTotNueBins[i])) << "%" << std::endl;
125  }
126  if(fTotNumu > 0) std::cout << "Numu efficiency = " << 100.*(fSelNumu/static_cast<float>(fTotNumu)) << "%" << std::endl;
127  for(unsigned int i = 0; i < fTotNumuBins.size(); ++i){
128  if(fTotNumuBins[i] > 0) std::cout << "Numu efficiency " << i << "= " << 100.*(fSelNumuBins[i]/static_cast<float>(fTotNumuBins[i])) << "%" << std::endl;
129  }
130  */
131  }
132 
133  //......................................................................
135  {
136 
137  /// Define containers for the things we're going to produce
138  std::unique_ptr< std::vector<Result> >
139  resultCol(new std::vector<Result>);
140 
141  /// Load in the pixel maps
142  std::vector< art::Ptr< cvn::PixelMap > > pixelmaplist;
144  auto pixelmapListHandle = evt.getHandle< std::vector< cvn::PixelMap > >(itag1);
145  if (pixelmapListHandle)
146  art::fill_ptr_vector(pixelmaplist, pixelmapListHandle);
147 
148  /// Make sure we have a valid name for the CVN type
149  /*
150  if(fCVNType == "Caffe"){
151  if(pixelmaplist.size()>0){
152  std::pair<const float*, const float*> pairedoutput= fCaffeHandler.Predict(*pixelmaplist[0]);
153  const float* output=pairedoutput.first;
154  //const float* features=pairedoutput.second;
155 
156  resultCol->emplace_back(output, fNOutput);
157 
158  }
159  }*/
160  if(fCVNType == "TF" || fCVNType == "Tensorflow" || fCVNType == "TensorFlow"){
161  // If we have a pixel map then use the TF interface to give us a prediction
162  if(pixelmaplist.size() > 0){
163 
164  std::vector< std::vector<float> > networkOutput = fTFHandler.Predict(*pixelmaplist[0]);
165  // cvn::Result can now take a vector of floats and works out the number of outputs
166  resultCol->emplace_back(networkOutput);
167 
168  /*
169  for(auto const& resaux: (*resultCol))
170  {
171  std::cout << "Is antineutrino: " << resaux.GetIsAntineutrinoProbability() << std::endl;
172  std::cout << "CC Numu: " << resaux.GetNumuProbability() << std::endl;
173  std::cout << "CC Nue: " << resaux.GetNueProbability() << std::endl;
174  std::cout << "CC Nutau: " << resaux.GetNutauProbability() << std::endl;
175  std::cout << "NC: " << resaux.GetNCProbability() << std::endl;
176  std::cout << "CC QE: " << resaux.GetQEProbability() << std::endl;
177  std::cout << "CC Res: " << resaux.GetResProbability() << std::endl;
178  std::cout << "CC DIS: " << resaux.GetDISProbability() << std::endl;
179  std::cout << "CC Other: " << resaux.GetOtherProbability() << std::endl;
180  std::cout << "0 Protons: " << resaux.Get0protonsProbability() << std::endl;
181  std::cout << "1 Protons: " << resaux.Get1protonsProbability() << std::endl;
182  std::cout << "2 Protons: " << resaux.Get2protonsProbability() << std::endl;
183  std::cout << ">2 Protons: " << resaux.GetNprotonsProbability() << std::endl;
184  std::cout << "0 Pions: " << resaux.Get0pionsProbability() << std::endl;
185  std::cout << "1 Pions: " << resaux.Get1pionsProbability() << std::endl;
186  std::cout << "2 Pions: " << resaux.Get2pionsProbability() << std::endl;
187  std::cout << ">3 Pions: " << resaux.GetNpionsProbability() << std::endl;
188  std::cout << "0 Pizeros: " << resaux.Get0pizerosProbability() << std::endl;
189  std::cout << "1 Pizeros: " << resaux.Get1pizerosProbability() << std::endl;
190  std::cout << "2 Pizeros: " << resaux.Get2pizerosProbability() << std::endl;
191  std::cout << ">2 Pizeros: " << resaux.GetNpizerosProbability() << std::endl;
192  std::cout << "0 Neutrons: " << resaux.Get0neutronsProbability() << std::endl;
193  std::cout << "1 Neutrons: " << resaux.Get1neutronsProbability() << std::endl;
194  std::cout << "2 Neutrons: " << resaux.Get2neutronsProbability() << std::endl;
195  std::cout << ">2 Neutrons: " << resaux.GetNneutronsProbability() << std::endl << std::endl;
196 
197  std::cout << "Is antineutrino: " << resaux.PredictedIsAntineutrino() << std::endl;
198  std::cout << "Predicted flavour: " << resaux.PredictedFlavour() << std::endl;
199  std::cout << "Predicted interaction: " << resaux.PredictedInteraction() << std::endl;
200  std::cout << "Predicted protons: " << resaux.PredictedProtons() << std::endl;
201  std::cout << "Predicted pions: " << resaux.PredictedPions() << std::endl;
202  std::cout << "Predicted pizeros: " << resaux.PredictedPizeros() << std::endl;
203  std::cout << "Predicted neutrons: " << resaux.PredictedNeutrons() << std::endl;
204  }
205  */
206 
207  // Classify other pixel maps if they exist
208  if(fMultiplePMs){
209  for(unsigned int p = 1; p < pixelmaplist.size(); ++p){
210  std::vector< std::vector<float> > output = fTFHandler.Predict(*pixelmaplist[p]);
211  resultCol->emplace_back(output);
212  }
213  }
214 
215  }
216  }
217  else{
218  mf::LogError("CVNEvaluator::produce") << "CVN Type not in the allowed list: Tensorflow" << std::endl;
219  mf::LogError("CVNEvaluator::produce") << "Exiting without processing events" << std::endl;
220  return;
221  }
222 /*
223 // Truth level debug code
224 // mf::LogInfo("CVNEvaluator::produce") << " Predicted: " << (*resultCol)[0].PredictedInteractionType() << std::endl;
225 
226  // Leigh: temporary testing code for performance
227 
228  InteractionType interaction = kOther;
229 
230  // * monte carlo
231  std::vector<art::Ptr<simb::MCTruth> > mclist;
232  auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >("generator");
233  if (mctruthListHandle)
234  art::fill_ptr_vector(mclist, mctruthListHandle);
235 
236  //unsigned short nmclist= mclist.size();
237 
238  //std::cout<<"mctruth: "<<nmclist<<std::endl;
239 
240  art::Ptr<simb::MCTruth> truth = mclist[0];
241  simb::MCNeutrino truthN=truth->GetNeutrino();
242  //truth = mclist[0];
243 
244  interaction = GetInteractionType(truthN);
245 
246  // Get the interaction vertex from the end point of the neutrino. This is
247  // because the start point of the lepton doesn't make sense for taus as they
248  // are decayed by the generator and not GEANT
249  TVector3 vtx = truthN.Nu().EndPosition().Vect();
250  bool isFid = (fabs(vtx.X())<310 && fabs(vtx.Y())<550 && vtx.Z()>50 && vtx.Z()<1244);
251 
252 
253 // if(isFid && pixelmaplist.size() > 0){
254 
255  unsigned int correctedInt = static_cast<unsigned int>(interaction);
256  if(correctedInt == 13) correctedInt = 12;
257 
258  if(isFid){
259  std::cout << "This fiducial event is a true " << truthN.Nu().PdgCode() << " and is it CC? " << (truthN.CCNC()==0) << " (" << correctedInt << ")" << std::endl;
260  float nueProb = (*resultCol)[0].GetNueProbability();
261  float numuProb = (*resultCol)[0].GetNumuProbability();
262  float nutauProb = (*resultCol)[0].GetNutauProbability();
263  float ncProb = (*resultCol)[0].GetNCProbability();
264  std::cout << "Summed probabilities " << numuProb << ", " << nueProb << ", " << nutauProb << ", " << ncProb << std::endl;
265  }
266 
267  unsigned int predInt = static_cast<unsigned int>((*resultCol)[0].PredictedInteractionType());
268 
269  ++fTotal;
270 // std::cout << " Truth :: Predicted = " << correctedInt << " :: " << predInt << " (" << numuProb << ", " << nueProb << ", " << nutauProb << ", " << ncProb << ")" << std::endl;
271  if((correctedInt >= 0 && correctedInt <=3) && (predInt >= 0 && predInt <= 3)){
272  ++fCorrect;
273  }
274  if((correctedInt >= 4 && correctedInt <= 7) && (predInt >= 4 && predInt <= 7)){
275  ++fCorrect;
276  }
277  if((correctedInt >= 8 && correctedInt <=11) && (predInt >= 8 && predInt <= 11)){
278  ++fCorrect;
279  }
280  if(correctedInt == 12 && predInt == 12){
281  ++fCorrect;
282  }
283 
284  if(correctedInt == predInt){
285  ++fFullyCorrect;
286  }
287 
288  if(abs(truthN.Nu().PdgCode()) == 12 && truthN.CCNC()==0){
289  ++fTotNue;
290  if(truthN.Nu().E() < 2.0) ++fTotNueBins[0];
291  if(truthN.Nu().E() >= 2.0 && truthN.Nu().E() <= 6.0) ++fTotNueBins[1];
292  if(truthN.Nu().E() > 6.0) ++fTotNueBins[2];
293 
294  if(nueProb > 0.7){
295  ++fSelNue;
296  if(truthN.Nu().E() < 2.0) ++fSelNueBins[0];
297  if(truthN.Nu().E() >= 2.0 && truthN.Nu().E() <= 6.0) ++fSelNueBins[1];
298  if(truthN.Nu().E() > 6.0) ++fSelNueBins[2];
299  }
300  }
301  if(abs(truthN.Nu().PdgCode()) == 14 && truthN.CCNC()==0){
302  ++fTotNumu;
303  if(truthN.Nu().E() < 2.0) ++fTotNumuBins[0];
304  if(truthN.Nu().E() >= 2.0 && truthN.Nu().E() <= 6.0) ++fTotNumuBins[1];
305  if(truthN.Nu().E() > 6.0) ++fTotNumuBins[2];
306 
307  if(numuProb > 0.5){
308  ++fSelNumu;
309  if(truthN.Nu().E() < 2.0) ++fSelNumuBins[0];
310  if(truthN.Nu().E() >= 2.0 && truthN.Nu().E() <= 6.0) ++fSelNumuBins[1];
311  if(truthN.Nu().E() > 6.0) ++fSelNumuBins[2];
312  }
313  }
314 
315  }
316 */ // End of truth level debug code
317 
318  evt.put(std::move(resultCol), fResultLabel);
319 
320  }
321 
323 } // end namespace cvn
324 ////////////////////////////////////////////////////////////////////////
325 
326 
327 
328 
329 
330 
331 
std::vector< std::vector< float > > Predict(const PixelMap &pm)
Return prediction arrays for PixelMap.
bool fMultiplePMs
Number of outputs fron neural net.
CVNEvaluator(fhicl::ParameterSet const &pset)
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
Wrapper for caffe::Net which handles construction and prediction.
Definition: TFNetHandler.h:22
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
PixelMap for CVN.
Utility class for truth labels.
std::vector< unsigned int > fSelNumuBins
std::string fCVNType
Can use Caffe or Tensorflow.
std::vector< unsigned int > fTotNueBins
unsigned int fFullyCorrect
Result for CVN.
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
def move(depos, offset)
Definition: depos.py:107
p
Definition: test.py:223
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
std::vector< unsigned int > fTotNumuBins
TFNetHandler for CVN.
void produce(art::Event &evt)
cvn::TFNetHandler fTFHandler
std::string fPixelMapInput
Module label for input pixel maps.
std::vector< unsigned int > fSelNueBins
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
QTextStream & endl(QTextStream &s)