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