NueSelLID_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Brief please!
3 /// \author Jianming Bian - bianjm@physics.umn.edu
4 ///////////////////////////////////////////////////////////////////////
5 #include "ShowerLID/EventLID.h"
6 #include "ShowerLID/ShowerLID.h"
7 #include "RecoJMShower/EID.h"
9 
10 extern "C" {
11 #include <sys/types.h>
12 #include <sys/stat.h>
13 }
14 #include <cmath>
15 #include <iostream>
16 #include <fstream>
17 #include <vector>
18 
19 // ROOT includes
20 #include "TFile.h"
21 #include "TH1D.h"
22 #include "TLorentzVector.h"
23 #include "TMath.h"
24 #include "TMultiLayerPerceptron.h"
25 #include "TRandom.h"
26 #include "TTree.h"
27 
28 // NOvA includes
29 #include "Calibrator/Calibrator.h"
30 #include "DAQChannelMap/DAQChannelMapConstants.h"
31 #include "Geometry/Geometry.h"
33 #include "RecoBase/Cluster.h"
34 #include "RecoBase/FilterList.h"
35 #include "RecoBase/Prong.h"
36 #include "RecoBase/RecoHit.h"
37 #include "RecoBase/Shower.h"
38 #include "RecoBase/Vertex.h"
39 #include "Utilities/AssociationUtil.h"
40 
41 // ART includes
49 #include "fhiclcpp/ParameterSet.h"
50 
51 
52 namespace jmshower
53 {
54 
55  class NueSelLID: public art::EDProducer
56  {
57  public:
58  explicit NueSelLID(const fhicl::ParameterSet& pset);
59  ~NueSelLID();
60 
61  virtual void produce(art::Event& evt);
62 
63  virtual void reconfigure(const fhicl::ParameterSet& pset);
64 
65  virtual void beginRun(art::Run& run);
66  virtual void endJob();
67 
68  double EnergyEstimator(TVector3 beamDir, TLorentzVector elecP4, double etot, int mode);
69  double GetPointDoca(TVector3 vtx, TVector3 start, TVector3 stop);
70 
71  protected:
72 
73  std::string fLIDLabel; ///< Label for LID input
74  std::string fSlicerLabel; ///< Label for rb::Cluster input
75  std::string fVertexLabel; ///< Label for rb::Vertex input
76  std::string fLibPath; ///< Path to PID libraries
77  // Osc. weight histograms
78 
79  TTree* fEvent;
80  //True Nuetrino info
86  double m_evtTrueNuP4[4];
87  double m_evtTrueNuVtx[3];
88  int m_evtTruePdg[20];
89  int m_evtTrueMother[20];
90  double m_evtTrueEnergy[20];
91  //Reco info
92  int m_evtRun;
93  int m_evtSRun;
95  double m_evtNSlice;
96  double m_evtNSh;
97  double m_evt;
98  double m_evtANN;
101  double m_evtEtot;
102  double m_evtSigOscW;
103  double m_evtBkgOscW;
112  double m_evtSh1Start[3];
113  double m_evtSh1Stop[3];
114  double m_evtSh1Dir[3];
115  double m_evtVtx[3];
116 
127  double m_evtSh1EPiLLL;
129  double m_evtSh1VtxGeV;
131  double m_evtSh1ShE;
132  double m_evtSh1Gap;
133  double m_evtSh1ELLL;
134  double m_evtSh1ELLT;
137  double m_evtSh1GLLL;
139 
144 
145 
146  private:
147 
148  TH1D *ht_nue_rw;
151 
152  TMultiLayerPerceptron *mlp_evt;
153 
154 
155  double OscWeight(int pdg, int ccnc, double energy, bool issig) ;
156  bool LoadHistograms(int detid);
157  void InitializeNTuple();
158  };
159 } // namespace
160 
161 
162 namespace jmshower
163 {
164  int Ncount0 = 0, Ncount1 = 0, Ncount2 = 0;
165  int Nfill = 0, Nput = 0;
166  // mass for electron, photon, muon, pi0, NA, proton, neutron and charged
167  const double xmass[7] = {0.000511, 0, 0.105658, 0.1349766, 0.938272, 0.939565379, 0.139570};
168  //......................................................................
170  : fEvent(0)
171  {
172  reconfigure(pset);
173  produces< std::vector<rb::Vertex> >();
174  produces< std::vector<jmshower::EID> >();
175  produces< art::Assns<jmshower::EID, rb::Cluster> >();
176  }
177 
178  //......................................................................
180  {
181  }
182 
183  //......................................................................
185  {
186  fLIDLabel = pset.get< std::string >("LIDLabel" );
187  fSlicerLabel = pset.get< std::string >("SlicerLabel");
188  fVertexLabel = pset.get< std::string >("VertexLabel");
189  fLibPath = pset.get< std::string >("LibraryPath");
190  }
191 
193  {
194  std::string fnameOs(fLibPath);
195  switch (detid) {
197  fnameOs += "extrahist_fd.root";
198  break;
199  case novadaq::cnv::kNDOS:
200  fnameOs += "extrahist_ndos.root";
201  break;
203  fnameOs += "extrahist_ndos.root";
204  break;
205  default:
206  throw cet::exception("BAD DETID") << __FILE__ << ":" << __LINE__
207  << " Bad detector id = " << detid;
208  }
209 /*
210  std::string fFnameOs;
211  cet::search_path sp("FW_SEARCH_PATH");
212  sp.find_file(fnameOs, fFnameOs);
213  struct stat sb;
214  if (fFnameOs.empty() || stat(fFnameOs.c_str(), &sb)!=0){
215  // Failed to resolve the file name
216  mf::LogWarning("NueSelLID") << "Histogram -- Os " << fFnameOs << " not found.";
217  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
218  }
219 */
220  TFile hfileOs(fnameOs.c_str());
221  ht_nue_rw = (TH1D*)hfileOs.Get("ht_nue_rw");
222  ht_bkg_nue_rw = (TH1D*)hfileOs.Get("ht_bkg_nue_rw");
223  ht_bkg_numu_rw = (TH1D*)hfileOs.Get("ht_bkg_numu_rw");
224 
225  std::string fnameWeightsEvt(fLibPath);
226  switch (detid) {
228  fnameWeightsEvt += "weights_evt_fd.txt";
229  break;
230  case novadaq::cnv::kNDOS:
231  fnameWeightsEvt += "weights_evt_ndos.txt";
232  break;
234  fnameWeightsEvt += "weights_evt_ndos.txt";
235  break;
236  default:
237  throw cet::exception("BAD DETID") << __FILE__ << ":" << __LINE__
238  << " Bad detector id = " << detid;
239  }
240 
241 // std::cout<<"Trained Evt ANN "<<fnameWeightsEvt.c_str()<<" is loaded."<<std::endl;
242  double egLLL;
243  double egLLT;
244  double epi0LLL;
245  double epi0LLT;
246  double epLLL;
247  double epLLT;
248  double enLLL;
249  double enLLT;
250  double epiLLL;
251  double epiLLT;
252  double pi0mass;
253  double vtxgev;
254  int type;
255  TTree *trainTree_evt = new TTree("trainTree_evt","m_trainTree_evt");
256  trainTree_evt->Branch("egLLL", &egLLL, "egLLL/D");
257  trainTree_evt->Branch("egLLT", &egLLT, "egLLT/D");
258  trainTree_evt->Branch("epi0LLL", &epi0LLL, "epi0LLL/D");
259  trainTree_evt->Branch("epi0LLT", &epi0LLT, "epi0LLT/D");
260  trainTree_evt->Branch("epLLL", &epLLL, "epLLL/D");
261  trainTree_evt->Branch("epLLT", &epLLT, "epLLT/D");
262  trainTree_evt->Branch("enLLL", &enLLL, "enLLL/D");
263  trainTree_evt->Branch("enLLT", &enLLT, "enLLT/D");
264  trainTree_evt->Branch("epiLLL", &epiLLL, "epiLLL/D");
265  trainTree_evt->Branch("epiLLT", &epiLLT, "epiLLT/D");
266  trainTree_evt->Branch("pi0mass", &pi0mass, "pi0mass/D");
267  trainTree_evt->Branch("vtxgev", &vtxgev, "vtxgev/D");
268  trainTree_evt->Branch("type", &type, "type/I");
269  mlp_evt = new TMultiLayerPerceptron("egLLL,egLLT,epi0LLL,epi0LLT,epLLL,epLLT,enLLL,enLLT,epiLLL,epiLLT,pi0mass,vtxgev:18:12:6:type",trainTree_evt);
270  mlp_evt->LoadWeights(fnameWeightsEvt.c_str());
271 
272  return true;
273  }
274 
275 
276  //......................................................................
278  {
279  if(fEvent) return; // Don't do this twice
280 
284 
285  //True neutrino information
286  fEvent = tfs->make<TTree>("fEvent", "fEvent");
287  fEvent->Branch("evtTrueNuCCNC", &m_evtTrueNuCCNC);
288  fEvent->Branch("evtTrueNuMode", &m_evtTrueNuMode);
289  fEvent->Branch("evtTrueNuPdg", &m_evtTrueNuPdg);
290  fEvent->Branch("evtTrueNuEnergy", &m_evtTrueNuEnergy);
291  fEvent->Branch("evtTrueNuIsFid", &m_evtTrueNuIsFid);
292  fEvent->Branch("evtTrueNuP4[4]", m_evtTrueNuP4);
293  fEvent->Branch("evtTrueNuVtx[3]", m_evtTrueNuVtx);
294  fEvent->Branch("evtTrueEnergy[20]", m_evtTrueEnergy);
295  fEvent->Branch("evtTruePdg[20]", m_evtTruePdg);
296  fEvent->Branch("evtTrueMother[20]", m_evtTrueMother);
297 
298  //Reconstructed information
299  fEvent->Branch("evtRun", &m_evtRun);
300  fEvent->Branch("evtSRun", &m_evtSRun);
301  fEvent->Branch("evtNSlice", &m_evtNSlice);
302  fEvent->Branch("evtNSh", &m_evtNSh);
303  fEvent->Branch("evtANN", &m_evtANN);
304  fEvent->Branch("evtANNSh1", &m_evtSh1DedxANN);
305  fEvent->Branch("evtANNESh1", &m_evtSh1DedxANNE);
306  fEvent->Branch("evtEtot", &m_evtEtot);
307  fEvent->Branch("evtSigOscW", &m_evtSigOscW);
308  fEvent->Branch("evtBkgOscW", &m_evtBkgOscW);
309  fEvent->Branch("evtIsMuon", &m_evtIsMuon);
310  fEvent->Branch("evtIsMuonSh1", &m_evtIsMuonSh1);
311  fEvent->Branch("evtIsInvPhoton", &m_evtIsInvPhoton);
312  fEvent->Branch("evtDetEnergy", &m_evtDetEnergy);
313  fEvent->Branch("evtNuEnergy", &m_evtNuEnergy);
314  fEvent->Branch("evtSh1Energy", &m_evtSh1Energy);
315  fEvent->Branch("evtSh1ExclEnergy", &m_evtSh1ExclEnergy);
316  fEvent->Branch("evtSh1Length", &m_evtSh1Length);
317  fEvent->Branch("evtSh1Start[3]", m_evtSh1Start);
318  fEvent->Branch("evtSh1Stop[3]", m_evtSh1Stop);
319  fEvent->Branch("evtSh1Dir[3]", m_evtSh1Dir);
320 
321  fEvent->Branch("evtVtx[3]", m_evtVtx);
322  fEvent->Branch("evtSh1EGLLL",&m_evtSh1EGLLL);
323  fEvent->Branch("evtSh1EGLLT",&m_evtSh1EGLLT);
324  fEvent->Branch("evtSh1EMuLLL",&m_evtSh1EMuLLL);
325  fEvent->Branch("evtSh1EMuLLT",&m_evtSh1EMuLLT);
326  fEvent->Branch("evtSh1EPi0LLL",&m_evtSh1EPi0LLL);
327  fEvent->Branch("evtSh1EPi0LLT",&m_evtSh1EPi0LLT);
328  fEvent->Branch("evtSh1EPLLL",&m_evtSh1EPLLL);
329  fEvent->Branch("evtSh1EPLLT",&m_evtSh1EPLLT);
330  fEvent->Branch("evtSh1ENLLL",&m_evtSh1ENLLL);
331  fEvent->Branch("evtSh1ENLLT",&m_evtSh1ENLLT);
332  fEvent->Branch("evtSh1EPiLLL",&m_evtSh1EPiLLL);
333  fEvent->Branch("evtSh1EPiLLT",&m_evtSh1EPiLLT);
334  fEvent->Branch("evtSh1VtxGeV",&m_evtSh1VtxGeV);
335  fEvent->Branch("evtSh1Pi0Mass",&m_evtSh1Pi0Mass);
336  fEvent->Branch("evtSh1ShE",&m_evtSh1ShE);
337  fEvent->Branch("evtSh1Gap",&m_evtSh1Gap);
338  fEvent->Branch("evtSh1ELLL",&m_evtSh1ELLL);
339  fEvent->Branch("evtSh1ELLT",&m_evtSh1ELLT);
340  fEvent->Branch("evtSh1MuLLL",&m_evtSh1MuLLL);
341  fEvent->Branch("evtSh1MuLLT",&m_evtSh1MuLLT);
342  fEvent->Branch("evtSh1GLLL",&m_evtSh1GLLL);
343  fEvent->Branch("evtSh1InvGLLL",&m_evtSh1InvGLLL);
344 
345  fEvent->Branch("evtSh1DistToEdge", &m_evtSh1DistToEdge);
346  fEvent->Branch("evtSh1SliceEdge2GeV", &m_evtSh1SliceEdge2GeV);
347  fEvent->Branch("evtSh1SliceEdge5GeV", &m_evtSh1SliceEdge5GeV);
348  fEvent->Branch("evtSh1SliceEdge10GeV", &m_evtSh1SliceEdge10GeV);
349 
350  }
351 
353  {
354  //fEvent
355  m_evtTrueNuCCNC = 0;
356  m_evtTrueNuMode = 0;
357  m_evtTrueNuPdg = 0;
358  m_evtTrueNuEnergy = 0;
359  m_evtTrueNuIsFid = 0;
360  for(int i=0;i<4;i++){m_evtTrueNuP4[i] = 0;}
361  for(int i=0;i<3;i++){m_evtTrueNuVtx[i] = 0;}
362  for(int i=0;i<20;i++){m_evtTruePdg[i] = 0;}
363  for(int i=0;i<20;i++){m_evtTrueMother[i] = 0;}
364  for(int i=0;i<20;i++){m_evtTrueEnergy[i] = 0;}
365  m_evtRun = 0;
366  m_evtSRun = 0;
367  m_evtEvent = 0;
368  m_evtNSlice = 0;
369  m_evtNSh = 0;
370  m_evt = 0;
371  m_evtANN = -5;
372  m_evtSh1DedxANN = -5;
373  m_evtSh1DedxANNE = -5;
374  m_evtEtot = 0;
375  m_evtSigOscW = 0;
376  m_evtBkgOscW = 0;
377  m_evtIsMuon = 0;
378  m_evtIsInvPhoton = 0;
379  m_evtDetEnergy = 0;
380  m_evtTrueNuEnergy = 0;
381  m_evtNuEnergy = 0;
382  m_evtSh1Energy = 0;
383  m_evtSh1ExclEnergy = 0;
384  m_evtSh1Length = 0;
385 
386  for(int i=0;i<3;i++){
387  m_evtVtx[i]=0;
388  m_evtSh1Start[i]=0;
389  m_evtSh1Stop[i]=0;
390  m_evtSh1Dir[i]=0;
391  }
392 
393  m_evtSh1EGLLL = 0;
394  m_evtSh1EGLLT = 0;
395  m_evtSh1EMuLLL = 0;
396  m_evtSh1EMuLLT = 0;
397  m_evtSh1EPi0LLL = 0;
398  m_evtSh1EPi0LLT = 0;
399  m_evtSh1EPLLL = 0;
400  m_evtSh1EPLLT = 0;
401  m_evtSh1ENLLL = 0;
402  m_evtSh1ENLLT = 0;
403  m_evtSh1EPiLLL = 0;
404  m_evtSh1EPiLLT = 0;
405  m_evtSh1VtxGeV = 0;
406  m_evtSh1Pi0Mass = 0;
407  m_evtSh1ShE = 0;
408  m_evtSh1Gap = 0;
409  m_evtSh1ELLL = 0;
410  m_evtSh1ELLT = 0;
411  m_evtSh1MuLLL = 0;
412  m_evtSh1MuLLT = 0;
413  m_evtSh1GLLL = 0;
414  m_evtSh1InvGLLL = 0;
415  m_evtSh1DistToEdge = 0;
419  }
420 
421  //......................................................................
422  double NueSelLID::OscWeight(int pdg, int ccnc, double energy, bool issig)
423  {
424  // Obtain oscillation weight for swapping MC
425  double weight = 1.0;
426  if(fabs(pdg)==12&&ccnc==0&&issig==true){
427  weight = ht_nue_rw->Interpolate(energy);
428  }
429  if(fabs(pdg)==14&&ccnc==0&&issig==false) weight = 0.1787396954;
430  return weight;
431  }
432 
433  //............................
434  double NueSelLID::EnergyEstimator(TVector3 beamDir, TLorentzVector elecP4, double etot, int mode)
435  {
436  // Neutrino energy estimator
437  double Ee = elecP4.E();
438  double Ehad = etot - Ee;
439  double EhadCor = 0.282525 + 1.0766*Ehad;
440  if(mode!=0&&mode!=1) return etot;
441  // EM shower has been calibrated in RecoJMShower, so there is no correction for Ee.
442  // The overall correction is for the assymetric convoluton effect.
443  double Enue = (Ee + EhadCor)/0.98343;
444  return Enue;
445  }
446 
447  //.......................................................................
448  double NueSelLID::GetPointDoca(TVector3 vtx, TVector3 start, TVector3 stop)
449  {
450  double d = 9999;
451  double denorm = (start-stop).Mag();
452  if(denorm<0.00001){d = (vtx-start).Mag();}else{
453  d = ((vtx-start).Cross(vtx-stop)).Mag()/denorm;
454  TVector3 dir = stop-start;
455  }
456  return d;
457  }
458 
459 
460  //......................................................................
462  {
463  std::cout<<"Total number of events "<<Ncount0<<std::endl;
464  std::cout<<"Muon like events "<<Ncount1<<std::endl;
465  std::cout<<"Nue CC events "<<Ncount2<<std::endl;
466  std::cout<<"Nfill, Nput "<<Nfill<<" "<<Nput<<std::endl;
467  }
468 
469  //......................................................................
470 
472  {
473  //----------------------------------------------------
474  // Run, subrun and event information
475  //----------------------------------------------------
476  Ncount0++;
477  int run = evt.run();
478  int srun = evt.subRun();
479  int event = evt.id().event();
480  // std::cout<<"In NueSel event "<<event<<std::endl;
481 
482 
483  // Neutrino MC truth
484  TLorentzVector evtTrueNuP4(0,0,0,0);
485  TVector3 evtTrueNuVtx(0,0,0);
486  int evtTrueNuCCNC = -1;
487  int evtTrueNuMode = -1;
488  int evtTrueNuPdg = 0;
489  double evtTrueNuEnergy = 0;
490 
491  /*
492  if(!evt.isRealData()){
493  for(int i=0;i<1;i++){
494  art::Handle< std::vector<simb::MCTruth> > mc;
495  try {
496  evt.getByLabel("generator", mc);
497  if (mc->size()!=1) break;
498  }
499  catch (...) {
500  break;
501  }
502  if ((*mc)[0].NeutrinoSet()) {
503  const simb::MCNeutrino& gen = (*mc)[0].GetNeutrino();
504  const simb::MCParticle& nu = gen.Nu();
505  evtTrueNuP4 = nu.Momentum();
506  evtTrueNuVtx.SetXYZ(nu.EndX(),nu.EndY(),nu.EndZ());
507  evtTrueNuCCNC = gen.CCNC() ;
508  evtTrueNuMode = gen.Mode() ;
509  evtTrueNuPdg = nu.PdgCode() ;
510  evtTrueNuEnergy = nu.Momentum().E() ;
511  }
512  }
513  }
514  */
515 
516 
517  // Showers reconstruction
518 
521 
522 
523  // Event information with PID
524  std::unique_ptr<std::vector<jmshower::EID> > eidcol(new std::vector<jmshower::EID>);
525  std::unique_ptr< art::Assns<jmshower::EID, rb::Cluster> > assns(new art::Assns<jmshower::EID, rb::Cluster>);
526 
527 
528 // EID product for analysis.
530  evt.getByLabel(fSlicerLabel, slicecol);
531 
532 
533 
534  for(unsigned int sliceIdx = 0; sliceIdx < slicecol->size(); ++sliceIdx){//Reconstruct EID slice by slice
535  if(rb::IsFiltered(evt,slicecol,sliceIdx)) continue;
536  const rb::Cluster& slice = (*slicecol)[sliceIdx];
537  if(slice.IsNoise()) continue;
538 
539  double evtGapTNS = 1e8;
540  for(unsigned int iic = 0; iic < slicecol->size(); ++iic){
541  if(iic==sliceIdx)continue;
542  const rb::Cluster& slice1 = (*slicecol)[iic];
543  if(slice1.IsNoise()) continue;
544  // if(slice1.NCell()<20)continue;
545  if(evtGapTNS > fabs(slice1.MeanTNS()-slice.MeanTNS())){
546  evtGapTNS = fabs(slice1.MeanTNS()-slice.MeanTNS());
547  }
548  }
549 
550 
552  m_evtRun = run;
553  m_evtSRun = srun;
554  m_evtEvent = event;
555 
556  m_evtTrueNuCCNC = evtTrueNuCCNC;
557  m_evtTrueNuMode = evtTrueNuMode;
558  m_evtTrueNuPdg = evtTrueNuPdg;
559 
560 
561  m_evtTrueNuEnergy = evtTrueNuEnergy;
562  m_evtTrueNuP4[0] = evtTrueNuP4.Px();
563  m_evtTrueNuP4[1] = evtTrueNuP4.Py();
564  m_evtTrueNuP4[2] = evtTrueNuP4.Pz();
565  m_evtTrueNuP4[3] = evtTrueNuP4.E();
566  m_evtTrueNuVtx[0] = evtTrueNuVtx[0];
567  m_evtTrueNuVtx[1] = evtTrueNuVtx[1];
568  m_evtTrueNuVtx[2] = evtTrueNuVtx[2];
569 
570  const double wx = (slice.NXCell() > 0) ? slice.W(slice.XCell(0).get()) : 0;
571  const double wy = (slice.NYCell() > 0) ? slice.W(slice.YCell(0).get()) : 0;
572  //calculate slce reco energy
573  double slicegev = 0;
574  for(unsigned int i=0;i<slice.NCell();i++){
575  const art::Ptr<rb::CellHit>& chit = slice.Cell(i);
576  double w = wx;
577  if(geom->Plane(chit->Plane())->View()==geo::kY)w=wy;
578 
579  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, w));
580 
581  double gev=0;
582  if(rhit.IsCalibrated())gev = rhit.GeV();
583  slicegev+=gev;
584  }
585 
586  art::FindOneP<slid::EventLID> foLID(slicecol, evt, fLIDLabel); //one lid object per slice
587  art::FindManyP<slid::ShowerLID> fmShwLID(slicecol, evt, fLIDLabel);
588  art::FindManyP<rb::Vertex> fmv(slicecol, evt, fVertexLabel); // get vertices
589 
590  TVector3 evtVtx(0,0,0);
591  if(fmv.isValid()){
592  std::vector<art::Ptr<rb::Vertex> > v = fmv.at(sliceIdx);
593  for (unsigned int j=0; j<v.size(); ++j) {
594  //store vertex coordinates
595  evtVtx.SetX(v[j]->GetX());
596  evtVtx.SetY(v[j]->GetY());
597  evtVtx.SetZ(v[j]->GetZ());
598  break;
599  }
600  }
601 
602  // Get slice level PID from LIDBuilder
603  jmshower::EID eid;
604  std::vector< art::Ptr<slid::ShowerLID> > shwlids;
605  std::vector<art::Ptr<rb::Shower> > showercol;
606  shwlids.clear();
607  showercol.clear();
609 
610  if(foLID.isValid())
611  {
612  cet::maybe_ref<art::Ptr<slid::EventLID> const> rlid(foLID.at(sliceIdx));
613  elid = rlid.ref();
614  //if the pointer is valid
615  if (elid){
616  if(fmShwLID.isValid())
617  {
618  shwlids = fmShwLID.at(sliceIdx);
619  art::FindOneP<rb::Shower> fos(shwlids, evt, fLIDLabel);
620  for(unsigned int shwId = 0; shwId < shwlids.size(); ++shwId)
621  {
622  cet::maybe_ref<art::Ptr<rb::Shower> const> rshw(fos.at(shwId));
623  art::Ptr<rb::Shower> shw = rshw.ref();
624  showercol.push_back(shw);
625  }
626  }
627  }
628  }
629 
630 
631  int evtNSh = showercol.size();
632  // std::cout<<"nshowercol, nshowerlids "<<showercol.size()<<" "<<shwlids.size()<<std::endl;
633  if(evtNSh<1){ // No shower be recontructed
634  fEvent->Fill(); // Fill NTuple
635  eid.SetANN(-5.);
636  eid.SetANNE(-5.);
637  eid.SetVtx(evtTrueNuVtx);
638  eid.SetNuEnergy(slicegev*1.9);
639  eid.SetDetEnergy(0);
640  eid.SetEvtNCell(slice.NCell());
641  eid.SetEvtSumCosTheta(0);
642  eid.SetEvtSumP(0);
643  eid.SetEvtSumPt(0);
644  eid.SetEvtSumP0(0);
645  eid.SetEvtSumPt0(0);
646 
647  eid.SetEvtEtot(0);
648  eid.SetEvtMinX(99999);
649  eid.SetEvtMinY(99999);
650  eid.SetEvtMinZ(99999);
651  eid.SetEvtMaxX(-99999);
652  eid.SetEvtMaxY(-99999);
653  eid.SetEvtMaxZ(-99999);
654  eid.SetEvtGapTNS(evtGapTNS);
655  eid.SetEvtNCellToEdge(10000);
656 
657  eid.SetSh1Energy(0);
658  eid.SetSh1ExclEnergy(0);
659  eid.SetSh1TotalLength(0);
660  TVector3 v3null(0,0,0);
661  eid.SetSh1Start(v3null);
662  eid.SetSh1Stop(v3null);
663  eid.SetSh1Dir(v3null);
664  eid.SetSh1VtxDoca(0);
665  eid.SetSh1Gap(0);
666  eid.SetSh1NPlane(0);
667  eid.SetSh1XNPlane(0);
668  eid.SetSh1YNPlane(0);
669  eid.SetSh1NCell(0);
670  eid.SetSh1XNCell(0);
671  eid.SetSh1YNCell(0);
672  eid.SetSh1PID(-1);
673  eid.SetSh1NCellToEdge(10000);
674 
675  eid.SetSh1Sh2Dang(0);
676  eid.SetSh2Energy(0);
677  eid.SetSh2ExclEnergy(0);
678  eid.SetSh2TotalLength(0);
679  eid.SetSh2Start(v3null);
680  eid.SetSh2Stop(v3null);
681  eid.SetSh2Dir(v3null);
682  eid.SetSh2VtxDoca(0);
683  eid.SetSh2Gap(0);
684  eid.SetSh2NPlane(0);
685  eid.SetSh2XNPlane(0);
686  eid.SetSh2YNPlane(0);
687  eid.SetSh2NCell(0);
688  eid.SetSh2XNCell(0);
689  eid.SetSh2YNCell(0);
690  eid.SetSh2PID(-1);
691  eid.SetSh2NCellToEdge(10000);
692 
693 
694  eid.SetIsMuon(0);
695  eid.SetIsInvPhoton(0);
696  eid.SetSigOscW(0);
697  eid.SetBkgOscW(0);
698 
699  eid.SetEGLLL(-10);
700  eid.SetEGLLT(-10);
701  eid.SetEMuLLL(-10);
702  eid.SetEMuLLT(-10);
703  eid.SetEPi0LLL(-10);
704  eid.SetEPi0LLT(-10);
705  eid.SetEPLLL(-10);
706  eid.SetEPLLT(-10);
707  eid.SetENLLL(-10);
708  eid.SetENLLT(-10);
709  eid.SetEPiLLL(-10);
710  eid.SetEPiLLT(-10);
711  eid.SetVtxGeV(0);
712  eid.SetPi0Mass(0);
713  eid.SetShE(0);
714  eid.SetGap(0);
715  eid.SetELLL(-10);
716  eid.SetELLT(-10);
717  eid.SetMuLLL(-10);
718  eid.SetMuLLT(-10);
719  eid.SetGLLL(-10);
720  eid.SetInvGLLL(-10);
721 
722  eid.SetSh1DistToEdge(0);
723  eid.SetEEdge2Cell(0);
724  eid.SetEEdge5Cell(0);
725  eid.SetEEdge10Cell(0);
726 
727 
728  eidcol->push_back(eid);
729  Nfill++;
730  fEvent->Fill(); // Fill NTuple
731  }
732  else{
733 
734  std::vector<rb::Prong> vXCluster;
735  std::vector<rb::Prong> vYCluster;
736  vXCluster.clear();
737  vYCluster.clear();
738 
739 
740  double evtMinX = 99999;
741  double evtMinY = 99999;
742  double evtMinZ = 99999;
743  double evtMaxX = -99999;
744  double evtMaxY = -99999;
745  double evtMaxZ = -99999;
746  for(unsigned int i=0;i<showercol.size();i++){
747  if((showercol[i]->Start())[0]<evtMinX)evtMinX = (showercol[i]->Start())[0];
748  if((showercol[i]->Stop())[0]<evtMinX)evtMinX = (showercol[i]->Stop())[0];
749  if((showercol[i]->Start())[1]<evtMinY)evtMinY = (showercol[i]->Start())[1];
750  if((showercol[i]->Stop())[1]<evtMinY)evtMinY = (showercol[i]->Stop())[1];
751  if((showercol[i]->Start())[2]<evtMinZ)evtMinZ = (showercol[i]->Start())[2];
752  if((showercol[i]->Stop())[2]<evtMinZ)evtMinZ = (showercol[i]->Stop())[2];
753 
754  if((showercol[i]->Start())[0]>evtMaxX)evtMaxX = (showercol[i]->Start())[0];
755  if((showercol[i]->Stop())[0]>evtMaxX)evtMaxX = (showercol[i]->Stop())[0];
756  if((showercol[i]->Start())[1]>evtMaxY)evtMaxY = (showercol[i]->Start())[1];
757  if((showercol[i]->Stop())[1]>evtMaxY)evtMaxY = (showercol[i]->Stop())[1];
758  if((showercol[i]->Start())[2]>evtMaxZ)evtMaxZ = (showercol[i]->Start())[2];
759  if((showercol[i]->Stop())[2]>evtMaxZ)evtMaxZ = (showercol[i]->Stop())[2];
760 
761  const art::PtrVector<rb::CellHit>& xcells = showercol[i]->XCells();
762  const art::PtrVector<rb::CellHit>& ycells = showercol[i]->YCells();
763  rb::Prong ClustX(xcells,showercol[i]->Start(),showercol[i]->Dir());
764  rb::Prong ClustY(ycells,showercol[i]->Start(),showercol[i]->Dir());
765  vXCluster.push_back(ClustX);
766  vYCluster.push_back(ClustY);
767  }
768 
769 
770  std::vector<TLorentzVector> showerpGam;// Set shower momentum with photon assumption
771  showerpGam.clear();
772  std::vector<double> vtxdoca;
773 
774  for(unsigned int i = 0; i < showercol.size(); i++){
775  double showereraw = shwlids[i]->ShowerEnergy();
776  TLorentzVector showerptrk(showercol[i]->Dir()[0]*showereraw,showercol[i]->Dir()[1]*showereraw,showercol[i]->Dir()[2]*showereraw,showereraw);
777  double vtxDoca = NueSelLID::GetPointDoca(evtVtx, showercol[i]->Start(), showercol[i]->Stop());
778  vtxdoca.push_back(vtxDoca);
779  showerpGam.push_back(showerptrk);
780  }
781  // Looping over showers
782  int ish1 = 0;
783  int ish1e = 0;
784  double sh1ann = -1;
785  double sh1eann = -1;
786  double sh1eanne = -1;
787  double sh1gev = 0;
788  double sh1egev = 0;
789  double sh1elength = 0;
790  double evtetot = 0;
791  std::vector<TLorentzVector> showerP4Col;
792  std::vector<TLorentzVector> showerP40Col;
793  std::vector<int> showerpid;
794  showerP4Col.clear();
795  showerP40Col.clear();
796  showerpid.clear();
797  double sh1egLLL = -10;
798  double sh1egLLT = -10;
799  double sh1emuLLL = -10;
800  double sh1emuLLT = -10;
801  double sh1epi0LLL = -10;
802  double sh1epi0LLT = -10;
803  double sh1epLLL = -10;
804  double sh1epLLT = -10;
805  double sh1enLLL = -10;
806  double sh1enLLT = -10;
807  double sh1epiLLL = -10;
808  double sh1epiLLT = -10;
809  double sh1vtxgev = 0;
810  double sh1pi0mass = 0;
811  double sh1shE = 0;
812  double sh1gap = 0;
813  double sh1eLLL = -10;
814  double sh1eLLT = -10;
815  double sh1muLLL = -10;
816  double sh1muLLT = -10;
817 
818 // int evtNCellToEdge = 10000;
819 
820  for(unsigned int i = 0; i < showercol.size(); i++){
821  // if(evtNCellToEdge>showercol[i]->NCellToEdge())evtNCellToEdge=showercol[i]->NCellToEdge();
822  double showerenergy = shwlids[i]->ShowerEnergy();
823  double showerlength = showercol[i]->TotalLength();
824  // ANN deterimantion
825  double showerann = elid->Value();
826  double showeranne = elid->ValueWithE();
827 
828  double egLLL = shwlids[i]->EGLLL();
829  double egLLT = shwlids[i]->EGLLT();
830  double emuLLL = shwlids[i]->EMuLLL();
831  double emuLLT = shwlids[i]->EMuLLT();
832  double epi0LLL = shwlids[i]->EPi0LLL();
833  double epi0LLT = shwlids[i]->EPi0LLT();
834  double epLLL = shwlids[i]->EPLLL();
835  double epLLT = shwlids[i]->EPLLT();
836  double enLLL = shwlids[i]->ENLLL();
837  double enLLT = shwlids[i]->EGLLT();
838  double epiLLL = shwlids[i]->EPiLLL();
839  double epiLLT = shwlids[i]->EPiLLT();
840  double vtxgev = shwlids[i]->VertexEnergy();
841  double pi0mass = shwlids[i]->Pi0mass();
842  double shE = shwlids[i]->ShowerEFrac();
843  double gap = shwlids[i]->Gap();
844 
845  std::map<int, float> tPartLongLL = shwlids[i]->fPartLongLL;
846  std::map<int, float> tPartTransLL = shwlids[i]->fPartTransLL;
847  std::map<int, float>::iterator itr = tPartLongLL.find(slid::DedxParticleType::kELECTRON);
848  double eLLL, eLLT, muLLL, muLLT;
849  if (itr != tPartLongLL.end()) eLLL = itr->second;
850  else eLLL = -5;
851  itr = tPartTransLL.find(slid::DedxParticleType::kELECTRON);
852  if (itr != tPartTransLL.end()) eLLT = itr->second;
853  else eLLT = -5;
854  itr = tPartLongLL.find(slid::DedxParticleType::kMUON);
855  if (itr != tPartLongLL.end()) muLLL = itr->second;
856  else muLLL = -5;
857  itr = tPartTransLL.find(slid::DedxParticleType::kMUON);
858  if (itr != tPartTransLL.end()) muLLT = itr->second;
859  else muLLT = -5;
860 
861  if(showerann > sh1ann){// Looking for the the shower with the largest electron ANN then set it's ANN as event ANN
862  ish1 = i;
863  sh1ann = showerann;
864  sh1gev = showerenergy;
865  }
866 
867  if(showerenergy > sh1egev){// Looking for the leading shower then set it's ANN as event ANN
868  ish1e = i;
869  sh1eann = showerann;
870  sh1eanne = showeranne;
871  sh1egev = showerenergy;
872  sh1elength = showerlength;
873 
874  sh1egLLL = egLLL;
875  sh1egLLT = egLLT;
876  sh1emuLLL = emuLLL;
877  sh1emuLLT = emuLLT;
878  sh1epi0LLL = epi0LLL;
879  sh1epi0LLT = epi0LLT;
880  sh1epLLL = epLLL;
881  sh1epLLT = epLLT;
882  sh1enLLL = enLLL;
883  sh1enLLT = enLLT;
884  sh1epiLLL = epiLLL;
885  sh1epiLLT = epiLLT;
886  sh1vtxgev = vtxgev;
887  sh1pi0mass = pi0mass;
888  sh1shE = shE;
889  sh1gap = gap;
890  sh1eLLL = eLLL;
891  sh1eLLT = eLLT;
892  sh1muLLL = muLLL;
893  sh1muLLT = muLLT;
894  }
895  evtetot+=showerenergy;
896 
897 //assign P4 based on shower PID, calculated by longitudinal+transverse likelihoods
898  TVector3 p3(showercol[i]->Dir()[0]*showerenergy,showercol[i]->Dir()[1]*showerenergy,showercol[i]->Dir()[2]*showerenergy);
899  double maxll = -9999;
900  int lltype = -1;
901  for( int itype = 0; itype <11; ++itype){
902  std::map<int, float> tempLLL;
903  std::map<int, float> tempLLT;
904  tempLLL.clear();
905  tempLLT.clear();
906  tempLLL = shwlids[i]->fPartLongLL;
907  tempLLT = shwlids[i]->fPartTransLL;
908  if(tempLLL[itype] + tempLLT[itype] > maxll){
909  lltype = itype;
910  maxll = tempLLL[itype] + tempLLT[itype];
911  }
912  }
913 
914  double m = xmass[0];
915  if(lltype>=0&&lltype<7) m=xmass[lltype];
916  double showerEtot0 = sqrt(p3.X()*p3.X()+p3.Y()*p3.Y()+p3.Z()*p3.Z());
917  double P3 = sqrt(showerEtot0*showerEtot0+2*showerEtot0*m);
918  double E = showerEtot0+m;
919 // std::cout<<"lltype, m, showere, showerE0, shower E "<<lltype<<" "<<maxll<<" "<<showerenergy<<" "<<showerEtot0<<" "<<showerEtot<<std::endl;
920  TLorentzVector showerP40(showercol[i]->Dir()[0]*showerEtot0,showercol[i]->Dir()[1]*showerEtot0,showercol[i]->Dir()[2]*showerEtot0, showerEtot0);
921  TLorentzVector showerP4(showercol[i]->Dir()[0]*P3, showercol[i]->Dir()[1]*P3, showercol[i]->Dir()[2]*P3, E);
922 //std::cout<<"showerP0, showerP "<<showerP40.Mag()<<" "<<showerP4.Mag()<<std::endl;
923 // TLorentzVector showerP4(showercol[i]->Dir()[0]*showerenergy, showercol[i]->Dir()[1]*showerenergy,showercol[i]->Dir()[2]*showerenergy,showerenergy);
924  showerP40Col.push_back(showerP40); // 4-Momentum of showers without pid
925  showerP4Col.push_back(showerP4); // 4-Momentum of showers
926  showerpid.push_back(lltype);// PID of showers
927  }
928 
929 // looking for the secondary shower
930 
931  int ish2 = -1;
932  double sh2gev = 0;
933 
934  for(unsigned int i = 0; i < showercol.size(); i++){
935  if(i==(unsigned int)ish1e)continue;
936  if(shwlids[i]->ShowerEnergy() > sh2gev){
937  ish2 = i;
938  sh2gev = shwlids[i]->ShowerEnergy();
939  }
940  }
941 
942 // calculate total momentum, Pt and direction for the event
943  TLorentzVector evtSumP4(0,0,0,0);
944  TLorentzVector evtSumP40(0,0,0,0);
945 
946  for(unsigned int ii=0;ii<showerP4Col.size();ii++){
947  evtSumP4+=showerP4Col[ii];
948  }
949 
950  for(unsigned int ii=0;ii<showerP40Col.size();ii++){
951  evtSumP40+=showerP40Col[ii];
952  }
953 
954  TVector3 nuDir = geom->NuMIBeamDirection(); //(-0.845e-3, 52.552e-3, 0.998618);//Beam dir
955  double evtSumTheta = evtSumP4.Angle(nuDir);
956  double evtSumPt = evtSumP4.Vect().Mag()*sin(evtSumTheta);
957  double evtSumP = evtSumP4.Vect().Mag();
958 
959  double evtSumTheta0 = evtSumP40.Angle(nuDir);
960  double evtSumPt0 = evtSumP40.Vect().Mag()*sin(evtSumTheta0);
961  double evtSumP0 = evtSumP40.Vect().Mag();
962 
963 
964 
965  bool ismuon = shwlids[ish1]->IsMuon(); // Muon like shower?
966  bool ismuone = shwlids[ish1e]->IsMuon(); // Muon like shower?
967 
968  if(sh1ann>0.84&&ismuon==true)Ncount1++;
969 // if( ismuon == true ) return; // remove muon
970 
971  TVector3 beamDir = evtTrueNuP4.Vect(); // True Neutrino 4-momentum
972 
973 // Osc illation weigt for sigal
974  double sigw = OscWeight(evtTrueNuPdg,evtTrueNuCCNC, evtTrueNuEnergy, true); // Signal
975  double bkgw = OscWeight(evtTrueNuPdg,evtTrueNuCCNC, evtTrueNuEnergy, false); // Background
976 
977 // Neutrino energy reconstruction
978  double evtNuEnergy = shwlids[ish1e]->NueEnergy();
979 
980  if(sh1ann>0.84&&sh1gev>0&&ismuon==false)Ncount2++;
981 // std::cout<<"N NueCC "<<Ncount2<<std::endl;
982  Nfill++;
983  fEvent->Fill(); // Fill NTuple
984 
985  eid.SetANN(sh1eann);
986  eid.SetANNE(sh1eanne);
987  eid.SetVtx(evtTrueNuVtx);
988  eid.SetNuEnergy(evtNuEnergy);
989  eid.SetDetEnergy(evtetot);
990 
991  eid.SetEvtNCell(slice.NCell());
992 // eid.SetEvtSumCosTheta(cos(evtSumTheta));
993  eid.SetEvtSumCosTheta(shwlids[ish1e]->CosTheta());
994  eid.SetEvtSumP(evtSumP);
995  eid.SetEvtSumPt(evtSumPt);
996  eid.SetEvtSumP0(evtSumP0);
997  eid.SetEvtSumPt0(evtSumPt0);
998 
999  eid.SetEvtEtot(evtetot);
1000  eid.SetEvtMinX(evtMinX);
1001  eid.SetEvtMinY(evtMinY);
1002  eid.SetEvtMinZ(evtMinZ);
1003  eid.SetEvtMaxX(evtMaxX);
1004  eid.SetEvtMaxY(evtMaxY);
1005  eid.SetEvtMaxZ(evtMaxZ);
1006  eid.SetEvtGapTNS(evtGapTNS);
1007 // eid.SetEvtNCellToEdge(evtNCellToEdge);
1008 
1009  eid.SetSh1Energy(sh1egev);
1010  eid.SetSh1ExclEnergy(shwlids[ish1e]->HadronicEnergy());
1011  eid.SetSh1TotalLength(sh1elength);
1012  eid.SetSh1Start(showercol[ish1e]->Start());
1013  eid.SetSh1Stop(showercol[ish1e]->Stop());
1014  eid.SetSh1Dir(showercol[ish1e]->Dir());
1015 
1016  eid.SetSh1VtxDoca(vtxdoca[ish1e]);
1017  eid.SetSh1Gap(shwlids[ish1e]->Gap());
1018  eid.SetSh1NPlane(showercol[ish1e]->ExtentPlane());
1019  eid.SetSh1XNPlane(vXCluster[ish1e].ExtentPlane());
1020  eid.SetSh1YNPlane(vYCluster[ish1e].ExtentPlane());
1021  eid.SetSh1NCell(showercol[ish1e]->NCell());
1022  eid.SetSh1XNCell(showercol[ish1e]->NXCell());
1023  eid.SetSh1YNCell(showercol[ish1e]->NYCell());
1024  eid.SetSh1PID(showerpid[ish1e]);
1025 // eid.SetSh1NCellToEdge(showercol[ish1e]->NCellToEdge());
1026  if(ish2!=-1){
1027  double evtSh1Sh2Dang = 180*showerP4Col[ish1e].Angle(showerP4Col[ish2].Vect())/3.14159265;
1028  eid.SetSh1Sh2Dang(evtSh1Sh2Dang);
1029  eid.SetSh2Energy(shwlids[ish2]->ShowerEnergy());
1030  eid.SetSh2ExclEnergy(shwlids[ish2]->HadronicEnergy());
1031  eid.SetSh2TotalLength(showercol[ish2]->TotalLength());
1032  eid.SetSh2Start(showercol[ish2]->Start());
1033  eid.SetSh2Stop(showercol[ish2]->Stop());
1034  eid.SetSh2Dir(showercol[ish2]->Dir());
1035 
1036  eid.SetSh2VtxDoca(vtxdoca[ish2]);
1037  eid.SetSh2Gap(shwlids[ish2]->Gap());
1038  eid.SetSh2NPlane(showercol[ish2]->ExtentPlane());
1039  eid.SetSh2XNPlane(vXCluster[ish2].ExtentPlane());
1040  eid.SetSh2YNPlane(vYCluster[ish2].ExtentPlane());
1041  eid.SetSh2NCell(showercol[ish2]->NCell());
1042  eid.SetSh2XNCell(showercol[ish2]->NXCell());
1043  eid.SetSh2YNCell(showercol[ish2]->NYCell());
1044  eid.SetSh2PID(showerpid[ish2]);
1045 // eid.SetSh2NCellToEdge(showercol[ish2]->NCellToEdge());
1046  }else{
1047  eid.SetSh1Sh2Dang(0);
1048  eid.SetSh2Energy(0);
1049  eid.SetSh2ExclEnergy(0);
1050  eid.SetSh2TotalLength(0);
1051  TVector3 v3null(0,0,0);
1052  eid.SetSh2Start(v3null);
1053  eid.SetSh2Stop(v3null);
1054  eid.SetSh2Dir(v3null);
1055  eid.SetSh2VtxDoca(0);
1056  eid.SetSh2Gap(0);
1057  eid.SetSh2NPlane(0);
1058  eid.SetSh2XNPlane(0);
1059  eid.SetSh2YNPlane(0);
1060  eid.SetSh2NCell(0);
1061  eid.SetSh2XNCell(0);
1062  eid.SetSh2YNCell(0);
1063  eid.SetSh2PID(-1);
1064  eid.SetSh2NCellToEdge(10000);
1065  }
1066 
1067 
1068  eid.SetIsMuon(ismuone);
1069  eid.SetIsInvPhoton(shwlids[ish1e]->IsInversePhoton());
1070  eid.SetSigOscW(sigw);
1071  eid.SetBkgOscW(bkgw);
1072 
1073  eid.SetEGLLL(sh1egLLL);
1074  eid.SetEGLLT(sh1egLLT);
1075  eid.SetEMuLLL(sh1emuLLL);
1076  eid.SetEMuLLT(sh1emuLLT);
1077  eid.SetEPi0LLL(sh1epi0LLL);
1078  eid.SetEPi0LLT(sh1epi0LLT);
1079  eid.SetEPLLL(sh1epLLL);
1080  eid.SetEPLLT(sh1epLLT);
1081  eid.SetENLLL(sh1enLLL);
1082  eid.SetENLLT(sh1enLLT);
1083  eid.SetEPiLLL(sh1epiLLL);
1084  eid.SetEPiLLT(sh1epiLLT);
1085  eid.SetVtxGeV(sh1vtxgev);
1086  eid.SetPi0Mass(sh1pi0mass);
1087  eid.SetShE(sh1shE);
1088  eid.SetGap(sh1gap);
1089  eid.SetELLL(sh1eLLL);
1090  eid.SetELLT(sh1eLLT);
1091  eid.SetMuLLL(sh1muLLL);
1092  eid.SetMuLLT(sh1muLLT);
1093 // eid.SetGLLL(shwlids[ish1e]->GLLL());
1094  eid.SetInvGLLL(shwlids[ish1e]->InverseGLLL());
1095 
1096 // eid.SetSh1DistToEdge(showercol[ish1e]->DistToEdge()); // Shower distance to edge
1097 // eid.SetEEdge2Cell(showercol[ish1e]->SliceEdge2GeV()); //Energy deposit within 2 cells to detector edges
1098 // eid.SetEEdge5Cell(showercol[ish1e]->SliceEdge5GeV()); //Energy deposit within 5 cells to detector edges
1099 // eid.SetEEdge10Cell(showercol[ish1e]->SliceEdge10GeV()); //Energy deposit within 10 cells to detector edges
1100  eidcol->push_back(eid);
1101  }
1102 
1103  util::CreateAssn(*this,evt,*eidcol,art::Ptr<rb::Cluster>(slicecol,sliceIdx),*assns);
1104  }
1105 
1106 
1107  //unsigned int iSlice = showercol[ish1e]->ISlice();
1108  //std::cout << "iSlice: " << iSlice << std::endl;
1109  //util::CreateAssn(*this,evt,*eidcol,slicelist[iSlice],*assns);
1110  //assert(iSlice >=0 && iSlice < slicecol->size());
1111 
1112 
1113  evt.put(std::move(eidcol));
1114  evt.put(std::move(assns));
1115  Nput++;
1116 // if(Ncount0%100==0)std::cout<<Ncount0<<" events processed, "<<Ncount2<<" Pass NueCC Sel "<<Ncount1<<" muon like"<<std::endl;
1117  return;
1118  }
1119 
1120 } //namespace jmshower
1121 
1122 
1123 ////////////////////////////////////////////////////////////////////////
1124 namespace jmshower
1125 {
1127 }
void SetEMuLLT(double emuLLT)
Definition: EID.cxx:663
void SetIsMuon(bool ismuon)
Definition: EID.cxx:543
SubRunNumber_t subRun() const
Definition: Event.h:72
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.
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Cluster.cxx:121
void SetSh1VtxDoca(double doca)
Definition: EID.cxx:299
void SetSh1YNPlane(unsigned int ynplane)
Definition: EID.cxx:334
void SetSh1NCellToEdge(int ncell)
Definition: EID.cxx:378
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void SetPi0Mass(double pi0mass)
Definition: EID.cxx:763
void SetEvtEtot(double etot)
Definition: EID.cxx:141
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
This is a helper class for ParticleIDAlg that provides a tidy structure in which to hold the dE/dx hi...
void SetEPiLLL(double epiLLL)
Definition: EID.cxx:733
void SetSh1Gap(double gap)
Definition: EID.cxx:307
void SetEPi0LLT(double epi0LLT)
Definition: EID.cxx:683
const Var weight
unsigned short Plane() const
Definition: CellHit.h:39
void SetSh2TotalLength(double sh1totallength)
Definition: EID.cxx:417
T sqrt(T number)
Definition: d0nt_math.hpp:156
bool LoadHistograms(int detid)
void SetSh2Stop(TVector3 sh1stop)
Definition: EID.cxx:436
void SetSh1XNPlane(unsigned int xnplane)
Definition: EID.cxx:325
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
void SetSh1NCell(unsigned int ncell)
Definition: EID.cxx:342
void SetSh1YNCell(unsigned int yncell)
Definition: EID.cxx:360
A collection of associated CellHits.
Definition: Cluster.h:47
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void SetNuEnergy(double energy)
Definition: EID.cxx:224
void SetSh2NCellToEdge(int ncell)
Definition: EID.cxx:534
const PlaneGeo * Plane(unsigned int i) const
void SetVtx(TVector3 vtx)
Definition: EID.cxx:581
void SetGap(double gap)
Definition: EID.cxx:783
DEFINE_ART_MODULE(TestTMapFile)
void SetGLLL(double gLLL)
Definition: EID.cxx:855
TODO.
Definition: FillPIDs.h:11
virtual void beginRun(art::Run &run)
double EnergyEstimator(TVector3 beamDir, TLorentzVector elecP4, double etot, int mode)
std::string fVertexLabel
Label for rb::Vertex input.
double Value() const
Definition: PID.h:22
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Definition: Run.h:31
::xsd::cxx::tree::type type
Definition: Database.h:110
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
void SetEPiLLT(double epiLLT)
Definition: EID.cxx:743
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void SetSh2XNPlane(unsigned int xnplane)
Definition: EID.cxx:481
void SetEvtSumP0(double p)
Definition: EID.cxx:103
void SetEPi0LLL(double epi0LLL)
Definition: EID.cxx:673
void SetEvtGapTNS(double tns)
Definition: EID.cxx:205
virtual void reconfigure(const fhicl::ParameterSet &pset)
Far Detector at Ash River, MN.
void SetSh2YNPlane(unsigned int ynplane)
Definition: EID.cxx:490
std::string fLIDLabel
Label for LID input.
void SetSh2YNCell(unsigned int yncell)
Definition: EID.cxx:516
void SetEvtMaxX(double m)
Definition: EID.cxx:178
void SetEEdge10Cell(double eedge10cell)
Definition: EID.cxx:904
void SetEvtSumP(double p)
Definition: EID.cxx:123
void SetSh2NPlane(unsigned int nplane)
Definition: EID.cxx:472
Float_t E
Definition: plot.C:20
void SetShE(double shE)
Definition: EID.cxx:773
void SetEvtMaxZ(double m)
Definition: EID.cxx:196
double GetPointDoca(TVector3 vtx, TVector3 start, TVector3 stop)
void SetSh2Start(TVector3 sh1start)
Definition: EID.cxx:427
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
void SetEvtMinY(double m)
Definition: EID.cxx:160
void SetSh2Dir(TVector3 sh1dir)
Definition: EID.cxx:445
Prototype Near Detector on the surface at FNAL.
void SetSh1NPlane(unsigned int nplane)
Definition: EID.cxx:316
T get(std::string const &key) const
Definition: ParameterSet.h:231
void SetELLL(double eLLL)
Definition: EID.cxx:793
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
void SetEGLLT(double egLLT)
Definition: EID.cxx:643
double energy
Definition: plottest35.C:25
void SetSh1PID(int pid)
Definition: EID.cxx:369
void SetVtxGeV(double vtxgev)
Definition: EID.cxx:753
float ValueWithE() const
Definition: EventLID.h:28
Near Detector in the NuMI cavern.
const double xmass[8]
Float_t d
Definition: plot.C:236
void SetSh2PID(int pid)
Definition: EID.cxx:525
void SetSh1XNCell(unsigned int xncell)
Definition: EID.cxx:351
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 SetSh1Sh2Dang(double sh1sh2dang)
Definition: EID.cxx:387
const double j
Definition: BetheBloch.cxx:29
double GetY(int dcm=1, int feb=0, int pix=0)
Definition: PlotOnMon.C:119
NueSelLID(const fhicl::ParameterSet &pset)
Vertex location in position and time.
std::string fSlicerLabel
Label for rb::Cluster input.
Definition: View.py:1
void SetEvtMaxY(double m)
Definition: EID.cxx:187
Definition: run.py:1
void SetSh2Energy(double sh1energy)
Definition: EID.cxx:398
OStream cout
Definition: OStream.cxx:6
void SetEPLLL(double epLLL)
Definition: EID.cxx:693
void SetEEdge2Cell(double eedge2cell)
Definition: EID.cxx:886
void SetEGLLL(double egLLL)
Definition: EID.cxx:633
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
EventNumber_t event() const
Definition: EventID.h:116
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
void SetEvtNCell(int ncell)
Definition: EID.cxx:84
void SetMuLLT(double muLLT)
Definition: EID.cxx:845
void SetSh1Dir(TVector3 sh1dir)
Definition: EID.cxx:289
float GeV() const
Definition: RecoHit.cxx:69
A Cluster with defined start position and direction.
Definition: Prong.h:19
T sin(T number)
Definition: d0nt_math.hpp:132
void SetInvGLLL(double invgLLL)
Definition: EID.cxx:865
double OscWeight(int pdg, int ccnc, double energy, bool issig)
T const * get() const
Definition: Ptr.h:321
TDirectory * dir
Definition: macro.C:5
void SetIsInvPhoton(bool isinvphoton)
Definition: EID.cxx:552
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void SetBkgOscW(double weight)
Definition: EID.cxx:571
void SetDetEnergy(double energy)
Definition: EID.cxx:75
void SetSh1ExclEnergy(double sh1exclenergy)
Definition: EID.cxx:252
void SetEvtSumCosTheta(double costheta)
Definition: EID.cxx:93
void SetEvtSumPt0(double pt)
Definition: EID.cxx:111
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
void SetSh1Start(TVector3 sh1start)
Definition: EID.cxx:271
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
void geom(int which=0)
Definition: geom.C:163
void SetMuLLL(double muLLL)
Definition: EID.cxx:835
virtual void produce(art::Event &evt)
float Mag() const
void SetEEdge5Cell(double eedge5cell)
Definition: EID.cxx:895
void SetSh2VtxDoca(double doca)
Definition: EID.cxx:455
void SetEvtMinZ(double m)
Definition: EID.cxx:169
void SetELLT(double eLLT)
Definition: EID.cxx:803
void SetSh1Energy(double sh1energy)
Definition: EID.cxx:242
void SetSh2XNCell(unsigned int xncell)
Definition: EID.cxx:507
std::string fLibPath
Path to PID libraries.
void SetENLLL(double enLLL)
Definition: EID.cxx:713
void SetEvtNCellToEdge(int ncell)
Definition: EID.cxx:215
void SetEvtMinX(double m)
Definition: EID.cxx:151
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
void SetANNE(double anne)
Definition: EID.cxx:21
void SetSh2ExclEnergy(double sh1exclenergy)
Definition: EID.cxx:408
RunNumber_t run() const
Definition: Event.h:77
Float_t w
Definition: plot.C:20
void SetSh1Stop(TVector3 sh1stop)
Definition: EID.cxx:280
void SetEMuLLL(double emuLLL)
Definition: EID.cxx:653
void SetEPLLT(double epLLT)
Definition: EID.cxx:703
void SetANN(double ann)
Definition: EID.cxx:12
EventID id() const
Definition: Event.h:56
void SetSh1DistToEdge(double dist)
Definition: EID.cxx:875
void SetENLLT(double enLLT)
Definition: EID.cxx:723
Encapsulate the geometry of one entire detector (near, far, ndos)
void SetEvtSumPt(double pt)
Definition: EID.cxx:132
double GetX(int ndb=14, int db=1, int feb=0, int pix=0)
Definition: PlotOnMon.C:111
TMultiLayerPerceptron * mlp_evt
TVector3 Vect() const
void SetSigOscW(double weight)
Definition: EID.cxx:561
void SetSh2NCell(unsigned int ncell)
Definition: EID.cxx:498
void SetSh2Gap(double gap)
Definition: EID.cxx:463
void SetSh1TotalLength(double sh1totallength)
Definition: EID.cxx:261