LIDTraining_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 // LIDTraining.cxx
3 //
4 // A module to reconstruct shower energy and topological variables
5 // for particle identification
6 //
7 // \author $Author: Jianming Bian $
8 ////////////////////////////////////////////////////////////////////////
9 #include <cmath>
10 #include <iostream>
11 #include <fstream>
12 #include <vector>
13 #include <algorithm>
14 
15 // ROOT includes
16 #include "TFile.h"
17 #include "TMath.h"
18 #include "TH1D.h"
19 #include "TNtuple.h"
20 #include "TFitter.h"
21 #include "TMinuit.h"
22 #include "TVector3.h"
23 #include "TGraph2DErrors.h"
24 
25 // NOvA includes
26 #include "Geometry/Geometry.h"
28 #include "GeometryObjects/Geo.h"
29 #include "Geometry/LiveGeometry.h"
30 #include "RecoBase/FilterList.h"
31 #include "RecoBase/Track.h"
32 #include "RecoBase/Vertex.h"
39 #include "Simulation/FLSHit.h"
40 
42 #include "Calibrator/Calibrator.h"
44 #include "MCCheater/BackTracker.h"
45 #include "RecoBase/RecoHit.h"
46 #include "RecoBase/Shower.h"
47 #include "ShowerLID/ShowerLID.h"
49 #include "ShowerLID/NuEEnergyAlg.h"
50 
51 // ART includes
63 
64 
65 #include "fhiclcpp/ParameterSet.h"
67 
68 #include "SummaryData/POTSum.h"
69 #include "TRandom.h"
70 #include "TMatrixD.h"
71 
72 #include "Utilities/func/MathUtil.h"
73 //#include "CosmicCalib/CalibUtil.h"
74 
75 
76 namespace slid{
78  {
79  public:
80  explicit LIDTraining(fhicl::ParameterSet const& pset);
81  ~LIDTraining();
82  virtual void analyze(const art::Event& evt);
83  virtual void beginJob();
84  virtual void beginSubRun(const art::SubRun& sr);
85  virtual void beginRun(const art::Run& run);
86  virtual void reconfigure(const fhicl::ParameterSet& pset);
87 
88  virtual void endJob();
89 
90  typedef std::vector<TLorentzVector> Vp4;
91  typedef std::vector<int> Vint;
92  typedef std::vector<double> Vdouble;
93 
94  TH1D *ht_nue_rw;
97 
98 //McTruth
99 
100  static const int GAMMAID;
101  static const int PIZEROID;
102  static const int PIPLUSID;
103  static const int PIMINUSID;
104  static const int NUEID;
105  static const int NUMUID;
106  static const int NUTAUID;
107 
108  private:
109  TTree* fEvent;
110  int m_evtRun;
114 
119  double m_evtTrueNuP4[4];
120  double m_evtTrueNuVtx[3];
121  double m_evtRwProb[20];
122  double m_evtRwRand[20];
124  double m_evtPi0Mgg;
127  double m_evtEtot;
128  double m_evtMinTop;
131  double m_evtMinBack;
132  double m_evtMinWest;
133  double m_evtMinEast;
134  double m_evtSh1Energy;
136  double m_evtSh1Length;
138  double m_evtSh1Gap;
139  double m_evtSh1IsFid;
141  double m_evtSh1P4[4];
142  double m_evtSh1Start[3];
143  double m_evtSh1Stop[3];
144  double m_evtSh1PlaneDedx[200];
146  double m_evtSh1DedxLLL[20];
147  double m_evtSh1DedxLLT[20];
152 
159  double m_evtSh1TrueP4[4];
174 
184 
187  double m_evtSh1SumP;
188  double m_evtSh2Energy;
189  double m_evtSheEnergy;
190 
191  TTree* fShower;
197  double m_showerTrueP4[4];
212 
213 
220 
221 
228  double m_showerStart[3];
229  double m_showerStop[3];
230  double m_showerDir[3];
231  double m_showerPlaneDedx[200];
237  double m_showerRun;
241  double m_showerId;
242  double m_showerEff;
243  double m_showerPur;
246 
248  unsigned int fMinTrack;
249  unsigned int fMaxTrack;
250  double fXThresh;
251  double fYThresh;
252  double fZMin;
253  double fZMax;
254  std::string fPotLabel; ///< Module that produced the POTSum object
257  bool fRecordMC;
266 
269  /// FCL parameters for NuE Energy alorithm
271  /// FCL parameters for particle ID alorithm (loglikelihoods and dE/dx)
273 
274 
275  };
276 }
277  //______________________________________________________________________________
278 namespace slid{
279  const int LIDTraining::GAMMAID = 22;
280  const int LIDTraining::PIZEROID = 111;
281  const int LIDTraining::PIPLUSID = 211;
282  const int LIDTraining::PIMINUSID = -211;
283  const int LIDTraining::NUEID = 12;
284  const int LIDTraining::NUMUID = 14;
285  const int LIDTraining::NUTAUID = 16;
286  const double xmass[7] = {0.000511, 0, 0.105658, 0.1349766, 0.938272, 0.939565379, 0.139570};
287  int Ncount0 = 0, Ncount1 = 0, Ncount2 = 0;
288 
290  EDAnalyzer(pset),
291  fParticleIDAlg(0),
292  fNuEEnergyAlg(0),
293  fNuEEnergyAlgPSet(pset.get< fhicl::ParameterSet >("NuEEnergyAlgPSet")),
294  fParticleIDAlgPSet(pset.get< fhicl::ParameterSet >("ParticleIDAlgPSet"))
295  {
296  std::cout<<"LIDTraining::LIDTraining(fhicl::ParameterSet const& pset)"<<std::endl;
297  reconfigure(pset);
298  }
299 
300 
301  //......................................................................
302 
304  {
305  fSliceLabel = pset.get< std::string >("SliceLabel" );
306  fProngLabel = pset.get< std::string >("ProngLabel" );
307  fInstLabel = pset.get< std::string >("InstLabel" );
308  fVertexLabel = pset.get< std::string >("VertexLabel" );
309  fHitLabel = pset.get< std::string >("HitLabel" );
310  fShowerLabel = pset.get< std::string >("ShowerLabel" );
311  fShowerLIDLabel = pset.get< std::string >("ShowerLIDLabel" );
312 
313  fCellHitPHThresh = pset.get< double >("CellHitPHThresh");
314  fMinTrack = pset.get< int >("MinTrack" );
315  fMaxTrack = pset.get< int >("MaxTrack" );
316  fXThresh = pset.get< double >("XThresh" );
317  fYThresh = pset.get< double >("YThresh" );
318  fZMin = pset.get< double >("ZMin" );
319  fZMax = pset.get< double >("ZMax" );
320  fPotLabel = pset.get< std::string >("PotLabel" );
321  fCutPlane = pset.get< int >("CutPlane" );
322  fRespectFilter = pset.get< bool >("RespectFilter" );
323  fTruthMatchByEnergy = pset.get< bool >("TruthMatchByEnergy");
324  fRecordMC = pset.get< bool >("RecordMC" );
325  fNuEEnergyAlgPSet = pset.get< fhicl::ParameterSet >("NuEEnergyAlgPSet");
326  fParticleIDAlgPSet = pset.get< fhicl::ParameterSet >("ParticleIDAlgPSet");
327  }
328 
329  //......................................................................
330 
332  {
333  if (fParticleIDAlg) delete fParticleIDAlg;
334  if (fNuEEnergyAlg) delete fNuEEnergyAlg;
335  }
336 
337  //......................................................................
338 
340  {
342  fEvent = tfs->make<TTree>("fEvent","fEvent");
343  //True neutrino information
344  fEvent->Branch("evtRun", &m_evtRun);
345  fEvent->Branch("evtSubRun", &m_evtSubRun);
346  fEvent->Branch("evtEvent", &m_evtEvent);
347  fEvent->Branch("evtSubEvent", &m_evtSubEvent);
348  fEvent->Branch("evtTrueNuCCNC", &m_evtTrueNuCCNC);
349  fEvent->Branch("evtTrueNuMode", &m_evtTrueNuMode);
350  fEvent->Branch("evtTrueNuPdg", &m_evtTrueNuPdg);
351  fEvent->Branch("evtTrueNuEnergy", &m_evtTrueNuEnergy);
352  fEvent->Branch("evtTrueNuP4[4]", m_evtTrueNuP4);
353  fEvent->Branch("evtTrueNuVtx[3]", m_evtTrueNuVtx);
354  fEvent->Branch("evtRwRand[20]", m_evtRwRand);
355  fEvent->Branch("evtRwProb[20]", m_evtRwProb);
356  fEvent->Branch("evtTotalRecoGeV", &m_evtTotalRecoGeV);
357  fEvent->Branch("evtPi0Mgg", &m_evtPi0Mgg);
358  fEvent->Branch("evtPi0StartDist", &m_evtPi0StartDist);
359  fEvent->Branch("evtPi0LineDist", &m_evtPi0LineDist);
360  fEvent->Branch("evtEtot", &m_evtEtot);
361  fEvent->Branch("evtMinTop", &m_evtMinTop);
362  fEvent->Branch("evtMinBottom", &m_evtMinBottom);
363  fEvent->Branch("evtMinFront", &m_evtMinFront);
364  fEvent->Branch("evtMinBack", &m_evtMinBack);
365  fEvent->Branch("evtMinEast", &m_evtMinEast);
366  fEvent->Branch("evtMinWest", &m_evtMinWest);
367 
368  fEvent->Branch("evtSh1Energy", &m_evtSh1Energy);
369  fEvent->Branch("evtNTrk", &m_evtNTrk);
370  fEvent->Branch("evtSh1DedxANN", &m_evtSh1DedxANN);
371  fEvent->Branch("evtSh1DedxANNEPi0", &m_evtSh1DedxANNEPi0);
372  fEvent->Branch("evtSh1DedxANNECos", &m_evtSh1DedxANNECos);
373  fEvent->Branch("evtSh1Length", &m_evtSh1Length);
374  fEvent->Branch("evtSh1NueEnergy", &m_evtSh1NueEnergy);
375  fEvent->Branch("evtSh1SumPt", &m_evtSh1SumPt);
376  fEvent->Branch("evtSh1SumP", &m_evtSh1SumP);
377  fEvent->Branch("evtSh1Gap", &m_evtSh1Gap);
378  fEvent->Branch("evtSh1IsFid", &m_evtSh1IsFid);
379  fEvent->Branch("evtSh1IsFidAna", &m_evtSh1IsFidAna);
380  fEvent->Branch("evtSh1P4[4]",m_evtSh1P4);
381  fEvent->Branch("evtSh1Start[3]",m_evtSh1Start);
382  fEvent->Branch("evtSh1Stop[3]",m_evtSh1Stop);
383  fEvent->Branch("evtSh1PlaneDedx[200]", m_evtSh1PlaneDedx);
384  fEvent->Branch("evtSh1CellPlaneDedx[66]", m_evtSh1CellPlaneDedx);
385  fEvent->Branch("evtSh1DedxLLL[20]", m_evtSh1DedxLLL);
386  fEvent->Branch("evtSh1DedxLLT[20]", m_evtSh1DedxLLT);
387  fEvent->Branch("evtSh1NPlane", &m_evtSh1NPlane);
388  fEvent->Branch("evtSh1NMIPPlane", &m_evtSh1NMIPPlane);
389  fEvent->Branch("evtSh1VtxGeV", &m_evtSh1VtxGeV);
390  fEvent->Branch("evtSh1Pi0Mgg", &m_evtSh1Pi0Mgg);
391  fEvent->Branch("evtSh1TrueDang", &m_evtSh1TrueDang);
392  fEvent->Branch("evtSh1TrueEDang", &m_evtSh1TrueEDang);
393  fEvent->Branch("evtSh1TrueP4[4]",m_evtSh1TrueP4);
394  fEvent->Branch("evtSh1TruePdg", &m_evtSh1TruePdg);
395  fEvent->Branch("evtSh1TruePur", &m_evtSh1TruePur);
396  fEvent->Branch("evtSh1TrueEff", &m_evtSh1TrueEff);
397  fEvent->Branch("evtSh1TrueDepEnergy", &m_evtSh1TrueDepEnergy);
398  fEvent->Branch("evtSh1TrueEEnergy", &m_evtSh1TrueEEnergy);
399  fEvent->Branch("evtSh1TrueGEnergy", &m_evtSh1TrueGEnergy);
400  fEvent->Branch("evtSh1TrueMuEnergy", &m_evtSh1TrueMuEnergy);
401  fEvent->Branch("evtSh1TruePEnergy", &m_evtSh1TruePEnergy);
402  fEvent->Branch("evtSh1TrueNEnergy", &m_evtSh1TrueNEnergy);
403  fEvent->Branch("evtSh1TruePiEnergy", &m_evtSh1TruePiEnergy);
404  fEvent->Branch("evtSh1TruePlaneEnergy[10]", m_evtSh1TruePlaneEnergy);
405  fEvent->Branch("evtSh1TrueEPlaneEnergy[10]", m_evtSh1TrueEPlaneEnergy);
406  fEvent->Branch("evtSh1TrueGPlaneEnergy[10]", m_evtSh1TrueGPlaneEnergy);
407  fEvent->Branch("evtSh1TrueMuPlaneEnergy[10]", m_evtSh1TrueMuPlaneEnergy);
408  fEvent->Branch("evtSh1TruePPlaneEnergy[10]", m_evtSh1TruePPlaneEnergy);
409  fEvent->Branch("evtSh1TrueNPlaneEnergy[10]", m_evtSh1TrueNPlaneEnergy);
410  fEvent->Branch("evtSh1TruePiPlaneEnergy[10]", m_evtSh1TruePiPlaneEnergy);
411 
412  fEvent->Branch("evtSh1TrueNavMotherPdg", &m_evtSh1TrueNavMotherPdg);
413  fEvent->Branch("evtSh1TrueNavMotherDang", &m_evtSh1TrueNavMotherDang);
414  fEvent->Branch("evtSh1TrueNavMotherId", &m_evtSh1TrueNavMotherId);
415  fEvent->Branch("evtSh1TrueNavMotherP4[4]", m_evtSh1TrueNavMotherP4);
416  fEvent->Branch("evtSh1TrueNavMotherPur", &m_evtSh1TrueNavMotherPur);
417  fEvent->Branch("evtSh1TrueNavMotherEff", &m_evtSh1TrueNavMotherEff);
418 
419  fEvent->Branch("evtSh1CosTheta", &m_evtSh1CosTheta);
420  fEvent->Branch("evtSh2Energy", &m_evtSh2Energy);
421 
422  fShower = tfs->make<TTree>("fShower","fShower");
423  fShower->Branch("showerId",&m_showerId);
424  fShower->Branch("showerRun",&m_showerRun);
425  fShower->Branch("showerSubRun",&m_showerSubRun);
426  fShower->Branch("showerEvent",&m_showerEvent);
427  fShower->Branch("showerSubEvent",&m_showerSubEvent);
428  fShower->Branch("showerNPlane", &m_showerNPlane);
429  fShower->Branch("showerStart[3]", m_showerStart);
430  fShower->Branch("showerStop[3]", m_showerStop);
431  fShower->Branch("showerDir[3]", m_showerDir);
432  fShower->Branch("showerEnergy", &m_showerEnergy);
433  fShower->Branch("showerNTrk", &m_showerNTrk);
434  fShower->Branch("showerIsSh1", &m_showerIsSh1);
435  fShower->Branch("showerDedxANN", &m_showerDedxANN);
436  fShower->Branch("showerDedxANNEPi0", &m_showerDedxANNEPi0);
437  fShower->Branch("showerDedxANNECos", &m_showerDedxANNECos);
438 
439  fShower->Branch("showerIsFid", &m_showerIsFid);
440  fShower->Branch("showerIsFidAna", &m_showerIsFidAna);
441  fShower->Branch("showerPlaneDedx[200]", m_showerPlaneDedx);
442  fShower->Branch("showerCellPlaneDedx[66]", m_showerCellPlaneDedx);
443  fShower->Branch("showerTransCellDedx[20]", m_showerTransCellDedx);
444 
445  fShower->Branch("showerTrueNuCCNC", &m_showerTrueNuCCNC);
446  fShower->Branch("showerTrueNuMode", &m_showerTrueNuMode);
447  fShower->Branch("showerTrueNuPdg", &m_showerTrueNuPdg);
448  fShower->Branch("showerTrueDang", &m_showerTrueDang);
449  fShower->Branch("showerTruePdg", &m_showerTruePdg);
450  fShower->Branch("showerTrueP4[4]", m_showerTrueP4);
451 
452  fShower->Branch("showerTrueDepEnergy", &m_showerTrueDepEnergy);
453  fShower->Branch("showerTrueEEnergy", &m_showerTrueEEnergy);
454  fShower->Branch("showerTrueGEnergy", &m_showerTrueGEnergy);
455  fShower->Branch("showerTrueMuEnergy", &m_showerTrueMuEnergy);
456  fShower->Branch("showerTruePEnergy", &m_showerTruePEnergy);
457  fShower->Branch("showerTrueNEnergy", &m_showerTrueNEnergy);
458  fShower->Branch("showerTruePiEnergy", &m_showerTruePiEnergy);
459  fShower->Branch("showerTruePlaneEnergy[10]", m_showerTruePlaneEnergy);
460  fShower->Branch("showerTrueEPlaneEnergy[10]", m_showerTrueEPlaneEnergy);
461  fShower->Branch("showerTrueGPlaneEnergy[10]", m_showerTrueGPlaneEnergy);
462  fShower->Branch("showerTrueMuPlaneEnergy[10]", m_showerTrueMuPlaneEnergy);
463  fShower->Branch("showerTruePPlaneEnergy[10]", m_showerTruePPlaneEnergy);
464  fShower->Branch("showerTrueNPlaneEnergy[10]", m_showerTrueNPlaneEnergy);
465  fShower->Branch("showerTruePiPlaneEnergy[10]", m_showerTruePiPlaneEnergy);
466 
467 
468 
469  fShower->Branch("showerEff", &m_showerEff);
470  fShower->Branch("showerPur", &m_showerPur);
471 
472  fShower->Branch("showerTrueNavMotherPdg", &m_showerTrueNavMotherPdg);
473  fShower->Branch("showerTrueNavMotherDang", &m_showerTrueNavMotherDang);
474  fShower->Branch("showerTrueNavMotherId", &m_showerTrueNavMotherId);
475  fShower->Branch("showerTrueNavMotherP4[4]", m_showerTrueNavMotherP4);
476  fShower->Branch("showerTrueNavMotherVtx[3]", m_showerTrueNavMotherVtx);
477  fShower->Branch("showerTrueNavMotherStop[3]", m_showerTrueNavMotherStop);
478  fShower->Branch("showerTrueNavmotherPur", &m_showerTrueNavMotherPur);
479  fShower->Branch("showerTrueNavmotherEff", &m_showerTrueNavMotherEff);
480 
481 
482  TFile hfileExtra("/nova/ana/users/bianjm/share/development/RecoJMShower/histograms/extrahist_fd.root");
483  ht_nue_rw = (TH1D*)hfileExtra.Get("ht_nue_rw");
484  ht_bkg_nue_rw = (TH1D*)hfileExtra.Get("ht_bkg_nue_rw");
485  ht_bkg_numu_rw = (TH1D*)hfileExtra.Get("ht_bkg_numu_rw");
486 
487  }
488 
490 /*
491  art::Handle< sumdata::POTSum > p;
492  sr.getByLabel(fPotLabel,p);
493  std::cout<<"pot = "<<p->totgoodpot<<std::endl;
494 */
495  }
496 
498  {
500  novadaq::cnv::DetId detId = fGeom->DetId();
501  if (fParticleIDAlg) delete fParticleIDAlg;
502  if (fNuEEnergyAlg) delete fNuEEnergyAlg;
505 
506  }
507 
508  //......................................................................
510  {
511 
512  Ncount0++;
513 
514  //----------------------------------------------------
515  // Run, subrun and event information
516  //----------------------------------------------------
517  int run = evt.run();
518  int srun = evt.subRun();
519  int event = evt.id().event();
520  std::cout<<"Event "<<event<<" ==============================================="<<std::endl;
522  evt.getByLabel(fSliceLabel, slicecol);
523  art::PtrVector<rb::Cluster> slicelist;
524  for(unsigned int i = 0; i < slicecol->size(); ++i){
525  slicelist.push_back(art::Ptr<rb::Cluster>(slicecol, i));
526  }
527 
529  const sim::ParticleNavigator & pnav = bt->ParticleNavigator();
530 
531  //Get all of the hits from calhits
533  evt.getByLabel(fHitLabel, hitsHandle);
535  for(unsigned int i = 0; i < hitsHandle->size(); ++i){
536  art::Ptr<rb::CellHit> hit(hitsHandle, i);
537  allHits.push_back(hit);
538  }
539  // We're gonna get a vector of NeutrinoEffPurs for allHits
540  std::vector<cheat::NeutrinoWithIndex> allNus = bt->allNeutrinoInteractions();
541 
542  std::vector< std::vector<cheat::NeutrinoEffPur> > sliceEffPurs;
543  sliceEffPurs.resize(slicelist.size());
544  for(unsigned int iSubvt = 0; iSubvt < slicelist.size(); ++iSubvt)
545  {
546  sliceEffPurs[iSubvt] = bt->SliceToMCTruth(slicelist[iSubvt]->AllCells(), allHits);//slicer4d
547 
548  }
549  std::vector<int> faveIds = bt->SliceToOrderedNuIdsByEff(allNus, sliceEffPurs);
550 
551  for(unsigned int sliceIdx = 0; sliceIdx < slicecol->size(); ++sliceIdx){//Reconstruct EID slice by slice
552  if(rb::IsFiltered(evt,slicecol,sliceIdx) && fRespectFilter) continue;
553  const rb::Cluster& slice = (*slicecol)[sliceIdx];
554  if(slice.IsNoise()) continue;
555 
556  if (fRecordMC && !(evt.isRealData()) &&(faveIds[sliceIdx] == -1)) continue;
557 
558  m_evtRun = run;
559  m_evtSubRun = srun;
560  m_evtEvent = event;
561  m_evtSubEvent = sliceIdx;
562  m_evtTrueNuEnergy = 0;
563  m_evtTrueNuP4[0] = 0;
564  m_evtTrueNuP4[1] = 0;
565  m_evtTrueNuP4[2] = 0;
566  m_evtTrueNuP4[3] = 0;
567  m_evtTrueNuVtx[0] = 0;
568  m_evtTrueNuVtx[1] = 0;
569  m_evtTrueNuVtx[2] = 0;
570 
571  for(int i=0;i<20;i++){
572  m_evtRwRand[i] = -1;
573  m_evtRwProb[i] = -1;
574  }
575 
576 
580 
581 
582  //-- ---------------------------------------------------------------------------------------------
583  // MC true information
584  //-----------------------------------------------------------------------------------------------
585  //Neutrino MC truth
586  TLorentzVector evtTrueNuP4(0,0,0,0);
587  TVector3 evtTrueNuVtx(0,0,-100);
588  int evtTrueNuCCNC = -1;
589  int evtTrueNuMode = -1;
590  int evtTrueNuPdg = 0;
591 
593  if(!evt.isRealData()&&fRecordMC) mc = sliceEffPurs[sliceIdx][0].neutrinoInt;
594  std::vector<const sim::Particle*> plist;
595  plist.clear();
596  if(!evt.isRealData()&&fRecordMC)plist = bt->MCTruthToParticles(mc);
597 
598  if(!evt.isRealData()&&fRecordMC){
599  if (mc->NeutrinoSet()) {
600  const simb::MCNeutrino& gen = mc->GetNeutrino();
601  const simb::MCParticle& nu = gen.Nu();
602  evtTrueNuP4 = nu.Momentum();
603  evtTrueNuVtx.SetXYZ(nu.EndX(),nu.EndY(),nu.EndZ());
604  evtTrueNuCCNC = gen.CCNC() ;
605  evtTrueNuMode = gen.Mode() ;
606  evtTrueNuPdg = nu.PdgCode() ;
607  }
608  }
609 
610 
611 
612  m_evtTrueNuCCNC = evtTrueNuCCNC;
613  m_evtTrueNuMode = evtTrueNuMode;
614  m_evtTrueNuPdg = evtTrueNuPdg;
615 
616 
617  m_evtTrueNuEnergy = evtTrueNuP4.E();
618  m_evtTrueNuP4[0] = evtTrueNuP4.Px();
619  m_evtTrueNuP4[1] = evtTrueNuP4.Py();
620  m_evtTrueNuP4[2] = evtTrueNuP4.Pz();
621  m_evtTrueNuP4[3] = evtTrueNuP4.E();
622  m_evtTrueNuVtx[0] = evtTrueNuVtx[0];
623  m_evtTrueNuVtx[1] = evtTrueNuVtx[1];
624  m_evtTrueNuVtx[2] = evtTrueNuVtx[2];
625 
626 
627  TRandom rand0(run*1100000000+srun*1000000+event);
628  double randnum0 = rand0.Uniform(0, 1);
629  TRandom rand1(run*1200000000+srun*100000+event);
630  double randnum1 = rand1.Uniform(0, 1);
631  TRandom rand2(run*1300000000+srun*1000000+event);
632  double randnum2 = rand2.Uniform(0, 1);
633 
634  double p_nue_rw = ht_nue_rw->Interpolate(evtTrueNuP4.E());
635  double p_bkg_nue_rw = 1;//ht_bkg_nue_rw->Interpolate(evtTrueNuP4.E());
636  double p_bkg_numu_rw = 1;//ht_bkg_numu_rw->Interpolate(evtTrueNuP4.E());
637 
638 
639  m_evtRwRand[0] = randnum0;
640  m_evtRwRand[1] = randnum1;
641  m_evtRwRand[2] = randnum2;
642  m_evtRwProb[0] = p_nue_rw;
643  m_evtRwProb[1] = p_bkg_nue_rw;
644  m_evtRwProb[2] = p_bkg_numu_rw;
645 
646  std::map<geo::OfflineChan, Vint> truePartPdg;
647  truePartPdg.clear();
648 
649  TLorentzVector electruep4(0,0,0,0);
650  std::vector<TLorentzVector> evtntpcol;
651  std::vector<TVector3> evtntvtxcol;
652  std::vector<TVector3> evtntpxzcol;
653  std::vector<TVector3> evtntpyzcol;
654  std::vector<int> evtntpdgcol;
655  std::vector<int> evtntmothercol;
656 
657  evtntpcol.clear();
658  evtntvtxcol.clear();
659  evtntpxzcol.clear();
660  evtntpyzcol.clear();
661  evtntpdgcol.clear();
662  evtntmothercol.clear();
663 
664  if(!evt.isRealData()&&fRecordMC){
665  for (unsigned int kk=0;kk<plist.size();kk++){
666  TLorentzVector ntp(0,0,0,0);
667  TVector3 ntvtx(0,0,0);
668  TVector3 ntpxz(0,0,0);
669  TVector3 ntpyz(0,0,0);
670  int ntpdg =0;
671  int ntmother = -1;
672  if(abs(plist[kk]->PdgCode()) == 11 && plist[kk]->Mother()==0){
673  electruep4 = plist[kk]->Momentum();
674  }
675  ntpdg = plist[kk]->PdgCode();
676  ntp = plist[kk]->Momentum();
677  ntpxz.SetXYZ(ntp.Px(),0,ntp.Pz());
678  ntpyz.SetXYZ(0,ntp.Py(),ntp.Pz());
679  ntmother = plist[kk]->Mother();
680  ntvtx.SetXYZ(plist[kk]->Vx(),plist[kk]->Vy(), plist[kk]->Vz());
681  evtntpcol.push_back(ntp);
682  evtntpxzcol.push_back(ntpxz);
683  evtntpyzcol.push_back(ntpyz);
684  evtntvtxcol.push_back(ntvtx);
685  evtntpdgcol.push_back(ntpdg);
686  evtntmothercol.push_back(ntmother);
687  }
688  }
689 
690  m_evtTotalRecoGeV = 0;
691  m_evtEtot = 0;
692  m_evtNTrk = 0;
693  m_evtMinTop = 9999;
694  m_evtMinBottom = 9999;
695 
696  m_evtMinFront = 9999;
697  m_evtMinBack = 9999;
698 
699  m_evtMinEast = 9999;
700  m_evtMinWest = 9999;
701 
702 
703  m_evtSh1DedxANN = 0;
706  m_evtPi0Mgg = 0;
707  m_evtPi0StartDist = 9999;
708  m_evtPi0LineDist = 9999;
709  m_evtSh1Energy = 0;
710  m_evtSh1Energy0 = 0;
711  m_evtSh1Length = 0;
712  m_evtSh1NueEnergy = 0;
713  m_evtSh1Gap = 0;
714  m_evtSh1P4[0] = 0;
715  m_evtSh1P4[1] = 0;
716  m_evtSh1P4[2] = 0;
717  m_evtSh1P4[3] = 0;
718  m_evtSh1Start[0] = 0;
719  m_evtSh1Start[1] = 0;
720  m_evtSh1Start[2] = 0;
721  m_evtSh1Stop[0] = 0;
722  m_evtSh1Stop[1] = 0;
723  m_evtSh1Stop[2] = 0;
724 
725  m_evtSh1TrueDang = 0;
726  m_evtSh1TrueEDang = 0;
727  m_evtSh1TruePdg = 0;
728  m_evtSh1TrueP4[0] = 0;
729  m_evtSh1TrueP4[1] = 0;
730  m_evtSh1TrueP4[2] = 0;
731  m_evtSh1TrueP4[3] = 0;
733  m_evtSh1TrueEEnergy = 0;
736  m_evtSh1TruePEnergy = 0;
737  m_evtSh1TrueNEnergy = 0;
739  for(unsigned int pl=0;pl<10;pl++){
747  }
748 
749 
750  m_evtSh1TruePur = 0;
751  m_evtSh1TrueEff = 0;
761 
762 
763  m_evtSh1IsFid = 0;
764  m_evtSh1IsFidAna = 0;
765  for(int itype=0;itype<20;itype++){
766  m_evtSh1DedxLLL[itype] = -9999;
767  m_evtSh1DedxLLT[itype] = -9999;
768  }
769 
770  for(int plane=0;plane<200;plane++){
772  }
773 
774  for(int cell=0;cell<11;cell++){
775  for(int plane=0;plane<6;plane++){
776  int cellplane = plane*11+cell;
777  m_evtSh1CellPlaneDedx[cellplane] = 0;
778  }
779  }
780  m_evtSh1NPlane = 0;
781  m_evtSh1NMIPPlane = 0;
782  m_evtSh1VtxGeV = 0;
783  m_evtSh1Pi0Mgg = 0;
784  m_evtSh1CosTheta = 0;
785  m_evtSh1SumPt=0;
786  m_evtSh1SumP=0;
787  m_evtSh2Energy = 0;
788  m_evtSheEnergy = 0;
789 
790  //-----------------------------------------------------
791  // Read slice, hit and track;
792  //-----------------------------------------------------
793 
794  std::map<geo::OfflineChan, int> hitslicemap;
795 
796  std::set<art::Ptr<rb::CellHit> > cellhitlist;
797  std::vector<const rb::Cluster* > cluster;
798 
799 
800  double totalRecoGeV = 0;
801  art::Ptr<rb::Cluster>clust(slicecol,sliceIdx);
802 
803  //get verticies associated with each slice
804  art::FindManyP<rb::Vertex> fmv(slicecol, evt, fVertexLabel);
805 
806  const double wx = (clust->NXCell() > 0) ? clust->W(clust->XCell(0).get()) : 0;
807  const double wy = (clust->NYCell() > 0) ? clust->W(clust->YCell(0).get()) : 0;
808  for(unsigned int i=0;i<clust->NXCell();i++){
809  const art::Ptr<rb::CellHit>& chit = clust->XCell(i);
810 
811  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, wx ));
812  if(rhit.IsCalibrated()){
813  if((clust->XCell(i)->PE())>fCellHitPHThresh ){
814  totalRecoGeV += rhit.GeV();
815  }
816  }
817  }
818  for(unsigned int i=0;i<clust->NYCell();i++){
819  const art::Ptr<rb::CellHit>& chit = clust->YCell(i);
820  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, wy));
821 
822  if(rhit.IsCalibrated()){
823  if((clust->YCell(i)->PE())>fCellHitPHThresh ){
824  totalRecoGeV += rhit.GeV();
825  }
826  }
827  }
828  m_evtTotalRecoGeV = totalRecoGeV;
829  //After read slice, hit and tracks.
830  // Determination of start point, end point and direction
831  // Fill true information of showers
832  std::vector<double> ntdangcol;
833  std::vector<TLorentzVector> ntpcol;
834  std::vector<TVector3> ntvtxcol;
835  std::vector<TVector3> ntstopcol;
836  std::vector<int> ntpdgcol;
837  std::vector<double> ntpurcol;
838  std::vector<double> nteffcol;
839 
840 
841  std::vector<TLorentzVector> ntnavmotherpcol;
842  std::vector<TVector3> ntnavmothervtxcol;
843  std::vector<TVector3> ntnavmotherstopcol;
844  std::vector<int> ntnavmotherpdgcol;
845  std::vector<double> ntnavmotherdangcol;
846  std::vector<int> ntnavmotheridcol;
847  std::vector<double> ntnavmotherpurcol;
848  std::vector<double> ntnavmothereffcol;
849 
850 
851 
852  ntdangcol.clear();
853  ntpcol.clear();
854  ntvtxcol.clear();
855  ntstopcol.clear();
856  ntpdgcol.clear();
857  ntpurcol.clear();
858  nteffcol.clear();
859 
860 
861  ntnavmotheridcol.clear();
862  ntnavmotherpdgcol.clear();
863  ntnavmotherdangcol.clear();
864  ntnavmotherpcol.clear();
865  ntnavmothervtxcol.clear();
866  ntnavmotherstopcol.clear();
867  ntnavmotherpurcol.clear();
868  ntnavmothereffcol.clear();
869 
870 
871 
872 
873  std::vector<art::Ptr<slid::ShowerLID> > slidcol;
874  slidcol.clear();
875  art::FindManyP<slid::ShowerLID> fmslid(slicecol, evt, fShowerLIDLabel);
876  if(fmslid.isValid()){
877  slidcol = fmslid.at(sliceIdx);
878  }
879  std::vector<const rb::Shower* > showercol;
880  showercol.clear();
881  art::FindMany<rb::Shower> fmshw(slicecol, evt, fShowerLabel);
882  if(fmshw.isValid()){
883  showercol = fmshw.at(sliceIdx);
884  }
885 
886  art::FindOneP<rb::Shower> fos(slidcol, evt, fShowerLIDLabel);
887 std::cout<<"run, srun, event, sevent, nsh =============================== "<<run<<" "<<srun<<" "<<event<<" "<<sliceIdx<<" "<<slidcol.size()<<std::endl;
888 
889 
890  m_evtNTrk = slidcol.size();
891 
892  int ish1 = -1;
893  int ish2 = -1;
894 
895  double sh1gev = 0;
896  double sh2gev = 0;
897 
898  for(int i = 0; i < (int)slidcol.size(); i++){
899  cet::maybe_ref<art::Ptr<rb::Shower> const> rshw(fos.at(i));
900  art::Ptr<rb::Shower> shw = rshw.ref();
901  double testE = fNuEEnergyAlg->ShowerEnergy(shw.get(), showercol, evt);
902  if(testE > sh1gev){
903  ish1 = i;
904  sh1gev = testE;
905  }
906  }
907 
908  for(int i = 0; i < (int)slidcol.size(); i++){
909  if(i==ish1)continue;
910  cet::maybe_ref<art::Ptr<rb::Shower> const> rshw(fos.at(i));
911  art::Ptr<rb::Shower> shw = rshw.ref();
912  double testE = fNuEEnergyAlg->ShowerEnergy(shw.get(), showercol, evt);
913  if(testE > sh2gev){
914  ish2 = i;
915  sh2gev = testE;
916  }
917  }
918 
919 
920  std::vector<double> showersumptcol;
921  std::vector<TLorentzVector> showersump4col;
922  showersumptcol.clear();
923  showersump4col.clear();
924 
925  for(unsigned int i = 0; i < slidcol.size(); i++){
926 
927  std::cout<<"Shower "<<i<<" ==============================================="<<std::endl;
928 
929  //For cell energy check and mc truth
930  cet::maybe_ref<art::Ptr<rb::Shower> const> rshw(fos.at(i));
931  art::Ptr<rb::Shower> shw = rshw.ref();
932  art::PtrVector< rb::CellHit > showerHits = shw->AllCells();
933  fParticleIDAlg->SetShower(shw.get(), showercol, evt);
934 
935  TVector3 trkdir=shw->Dir();
936  TVector3 trkstart=shw->Start();
937  TVector3 trkstop=shw->Stop();
938  double ntdang= 200.0;
939  TLorentzVector ntp(0,0,0,0);
940  TVector3 ntvtx(0,0,0);
941  TVector3 ntstop(0,0,0);
942  int ntpdg =0;
943 
944  int ntnavmotherpdg = 0;
945  double ntnavmotherdang = 0;
946  int ntnavmotherid = -1;
947  TLorentzVector ntnavmotherp(0,0,0,0);
948  TVector3 ntnavmothervtx(0,0,0);
949  TVector3 ntnavmotherstop(0,0,0);
950 
951  double ntpur = 0;
952  double nteff = 0;
953  double ntnavmotherpur = 0;
954  double ntnavmothereff = 0;
955 
956  int bestang =-1;
957  if(!evt.isRealData()&&fRecordMC){
958 
959  if( fTruthMatchByEnergy ){
960 
961  const std::vector< cheat::TrackIDE > trackides = bt->HitsToTrackIDE( showerHits );
962  if( !trackides.empty() ){
963  const sim::Particle * part = bt->TrackIDToParticle( trackides[0].trackID );
964  ntpdg = part->PdgCode();
965  double angle = part->Momentum().Vect().Angle( trkdir ) ;
966  double backangle = part->Momentum().Vect().Angle( -trkdir ) ;
967  ntdang = std::min( angle, backangle );
968  std::set< int > setid;
969  setid.insert( trackides[0].trackID );
970 
971  std::map<int, double>* purMap = new std::map<int, double>;
972  std::map<int, int>* parents = new std::map<int, int>;
973  ntpur = bt->HitCollectionPurity( setid, showerHits, purMap, parents, true);
974  nteff = bt->HitCollectionEfficiency( setid, showerHits, slicelist[sliceIdx]->AllCells(), geo::kXorY);
975  purMap->clear();
976  parents->clear();
977  delete purMap;
978  delete parents;
979 
980  m_showerPur = ntpur;
981  m_showerEff = nteff;
982  std::cout<<"Pdg code from energy "<<ntpdg<<std::endl;
983  }
984  else ntpdg = -1;
985  }// end if truthMatchByEnergy
986  else{
987  int partid;
988  for (unsigned int kk=0;kk<plist.size();kk++){
989  int abspdg = abs(plist[kk]->PdgCode());
990  if(abspdg==12||abspdg==14||abspdg==16)continue;
991 
992  TLorentzVector true_4mom = plist[kk]->Momentum();
993  if(true_4mom.Vect().Mag()<0.000001)continue;
994  double ntangd = trkdir.Angle(true_4mom.Vect());
995  ntangd=fmod(ntangd, TMath::Pi());
996  // if(ntangd>TMath::Pi()/2)
997  // ntangd = TMath::Pi()-ntangd;
998  if(ntangd<ntdang){
999  ntdang = ntangd;
1000  bestang = plist[kk]->PdgCode();
1001  ntpdg = plist[kk]->PdgCode();
1002  ntp = true_4mom;
1003  ntvtx.SetXYZ(plist[kk]->Vx(),plist[kk]->Vy(),plist[kk]->Vz());
1004  ntstop.SetXYZ(plist[kk]->EndX(),plist[kk]->EndY(),plist[kk]->EndZ());
1005  partid = plist[kk]->TrackId();
1006 
1007  int navmotherid = pnav.EveId(plist[kk]->TrackId());
1008  for (unsigned int ll=0;ll<plist.size();ll++){
1009  if(plist[ll]->TrackId()==navmotherid){
1010  ntnavmotherpdg = plist[ll]->PdgCode();
1011  ntnavmotherdang = trkdir.Angle(plist[ll]->Momentum().Vect());
1012  ntnavmotherdang = fmod(ntnavmotherdang, 2*3.14159265);
1013  ntnavmotherdang = ntnavmotherdang*180/3.14159265;
1014  ntnavmotherid = plist[ll]->TrackId();
1015  ntnavmotherp = plist[ll]->Momentum();
1016  ntnavmothervtx.SetXYZ(plist[ll]->Vx(),plist[ll]->Vy(),plist[ll]->Vz());
1017  ntnavmotherstop.SetXYZ(plist[ll]->EndX(),plist[ll]->EndY(),plist[ll]->EndZ());
1018  break;
1019  }
1020  }
1021  }
1022  }
1023  double showerTrueDepEnergy = 0;
1024  double showerTrueEEnergy = 0;
1025  double showerTrueGEnergy = 0;
1026  double showerTrueMuEnergy = 0;
1027  double showerTruePEnergy = 0;
1028  double showerTrueNEnergy = 0;
1029  double showerTruePiEnergy = 0;
1030  double showerTruePlaneEnergy[10];
1031  double showerTrueEPlaneEnergy[10];
1032  double showerTrueGPlaneEnergy[10];
1033  double showerTrueMuPlaneEnergy[10];
1034  double showerTruePPlaneEnergy[10];
1035  double showerTrueNPlaneEnergy[10];
1036  double showerTruePiPlaneEnergy[10];
1037 
1038  for(unsigned int pl=0;pl<10;pl++){
1039  showerTruePlaneEnergy[pl] = 0;
1040  showerTrueEPlaneEnergy[pl] = 0;
1041  showerTrueGPlaneEnergy[pl] = 0;
1042  showerTrueMuPlaneEnergy[pl] = 0;
1043  showerTruePPlaneEnergy[pl] = 0;
1044  showerTrueNPlaneEnergy[pl] = 0;
1045  showerTruePiPlaneEnergy[pl] = 0;
1046  }
1047  for(unsigned int ii=0;ii<showerHits.size();ii++){
1048  unsigned int cellplane = showerHits[ii]->Plane();
1049  std::vector<sim::FLSHit> tHits = bt->HitToFLSHit(showerHits[ii]);
1050 // std::cout<<"cellE nflsHits "<<showerHits[ii]->PE()<<" "<<tHits.size()<<std::endl;
1051  for(unsigned int n = 0; n < tHits.size(); n++){
1052  sim::FLSHit fls = tHits[n];
1053  showerTrueDepEnergy += fls.GetEdep();
1054  if(fabs(fls.GetPDG())==11){
1055  showerTrueEEnergy += fls.GetEdep();
1056  }
1057  if(fls.GetPDG()==22){
1058  showerTrueGEnergy += fls.GetEdep();
1059  }
1060  if(fabs(fls.GetPDG())==13){
1061  showerTrueMuEnergy += fls.GetEdep();
1062  }
1063  if(fabs(fls.GetPDG())==2212){
1064  showerTruePEnergy += fls.GetEdep();
1065  }
1066  if(fabs(fls.GetPDG())==2112){
1067  showerTrueNEnergy += fls.GetEdep();
1068  }
1069  if(fabs(fls.GetPDG())==211){
1070  showerTruePiEnergy += fls.GetEdep();
1071  }
1072 
1073 // for(unsigned int pl=0;pl<shw->ExtentPlane();pl++){
1074 // std::cout<<"cellplane , global plane "<<cellplane<<" "<<fParticleIDAlg->PlaneToGlobalPlane(pl)<<std::endl;
1075 // }
1076 
1077  for(unsigned int pl=0;pl<10;pl++){
1078  if(cellplane==fParticleIDAlg->PlaneToGlobalPlane(pl)){
1079  showerTruePlaneEnergy[pl]+=fls.GetEdep();
1080  if(fabs(fls.GetPDG())==11){
1081  showerTrueEPlaneEnergy[pl] += fls.GetEdep();
1082  }
1083  if(fls.GetPDG()==22){
1084  showerTrueGPlaneEnergy[pl] += fls.GetEdep();
1085  }
1086  if(fabs(fls.GetPDG())==13){
1087  showerTrueMuPlaneEnergy[pl] += fls.GetEdep();
1088  }
1089  if(fabs(fls.GetPDG())==2212){
1090  showerTruePPlaneEnergy[pl] += fls.GetEdep();
1091  }
1092  if(fabs(fls.GetPDG())==2112){
1093  showerTrueNPlaneEnergy[pl] += fls.GetEdep();
1094  }
1095  if(fabs(fls.GetPDG())==211){
1096  showerTruePiPlaneEnergy[pl] += fls.GetEdep();
1097  }
1098  }
1099  }
1100 
1101  }
1102  }
1103  m_showerTrueDepEnergy = showerTrueDepEnergy;
1104  m_showerTrueEEnergy = showerTrueEEnergy;
1105  m_showerTrueGEnergy = showerTrueGEnergy;
1106  m_showerTrueMuEnergy = showerTrueMuEnergy;
1107  m_showerTruePEnergy = showerTruePEnergy;
1108  m_showerTrueNEnergy = showerTrueNEnergy;
1109  m_showerTruePiEnergy = showerTruePiEnergy;
1110  for(unsigned int pl=0;pl<10;pl++){
1111  m_showerTruePlaneEnergy[pl] = showerTruePlaneEnergy[pl];
1112  m_showerTrueEPlaneEnergy[pl] = showerTrueEPlaneEnergy[pl];
1113  m_showerTrueGPlaneEnergy[pl] = showerTrueGPlaneEnergy[pl];
1114  m_showerTrueMuPlaneEnergy[pl] = showerTrueMuPlaneEnergy[pl];
1115  m_showerTruePPlaneEnergy[pl] = showerTruePPlaneEnergy[pl];
1116  m_showerTrueNPlaneEnergy[pl] = showerTrueNPlaneEnergy[pl];
1117  m_showerTruePiPlaneEnergy[pl] = showerTruePiPlaneEnergy[pl];
1118  }
1119 
1120  if((int)i==ish1){
1121  m_evtSh1TrueDepEnergy = showerTrueDepEnergy;
1122  m_evtSh1TrueEEnergy = showerTrueEEnergy;
1123  m_evtSh1TrueGEnergy = showerTrueGEnergy;
1124  m_evtSh1TrueMuEnergy = showerTrueMuEnergy;
1125  m_evtSh1TruePEnergy = showerTruePEnergy;
1126  m_evtSh1TrueNEnergy = showerTrueNEnergy;
1127  m_evtSh1TruePiEnergy = showerTruePiEnergy;
1128  for(unsigned int pl=0;pl<10;pl++){
1129  m_evtSh1TruePlaneEnergy[pl] = showerTruePlaneEnergy[pl];
1130  m_evtSh1TrueEPlaneEnergy[pl] = showerTrueEPlaneEnergy[pl];
1131  m_evtSh1TrueGPlaneEnergy[pl] = showerTrueGPlaneEnergy[pl];
1132  m_evtSh1TrueMuPlaneEnergy[pl] = showerTrueMuPlaneEnergy[pl];
1133  m_evtSh1TruePPlaneEnergy[pl] = showerTruePPlaneEnergy[pl];
1134  m_evtSh1TrueNPlaneEnergy[pl] = showerTrueNPlaneEnergy[pl];
1135  m_evtSh1TruePiPlaneEnergy[pl] = showerTruePiPlaneEnergy[pl];
1136  }
1137  }
1138 
1139  std::cout<<"Pdg code from angle "<<bestang<<std::endl;
1140  std::set< int > setid;
1141  setid.insert( partid );
1142  std::map<int, double>* purMap = new std::map<int, double>;
1143  std::map<int, int>* parents = new std::map<int, int>;
1144  ntpur = bt->HitCollectionPurity( setid, showerHits, purMap, parents, true);
1145  nteff = bt->HitCollectionEfficiency( setid, showerHits, slicelist[sliceIdx]->AllCells(), geo::kXorY);
1146  purMap->clear();
1147  parents->clear();
1148  delete purMap;
1149  delete parents;
1150  m_showerPur = ntpur;
1151  m_showerEff = nteff;
1152  std::set< int > setnavmotherid;
1153  setnavmotherid.insert( ntnavmotherid );
1154  std::map<int, double>* navmotherpurMap = new std::map<int, double>;
1155  std::map<int, int>* navmotherparents = new std::map<int, int>;
1156  ntnavmotherpur = bt->HitCollectionPurity( setnavmotherid, showerHits, navmotherpurMap, navmotherparents, true);
1157  ntnavmothereff = bt->HitCollectionEfficiency( setnavmotherid, showerHits, slicelist[sliceIdx]->AllCells(), geo::kXorY);
1158  navmotherpurMap->clear();
1159  navmotherparents->clear();
1160  delete navmotherpurMap;
1161  delete navmotherparents;
1162 
1163  m_showerTrueNavMotherPur = ntnavmotherpur;
1164  m_showerTrueNavMotherEff = ntnavmothereff;
1165  }
1166  }
1167  ntdang = ntdang*TMath::RadToDeg();
1168  std::cout<<"dangle "<<ntdang<<std::endl;
1169  ntdangcol.push_back(ntdang);
1170  ntpcol.push_back(ntp);
1171  ntvtxcol.push_back(ntvtx);
1172  ntstopcol.push_back(ntstop);
1173  ntpdgcol.push_back(ntpdg);
1174  ntpurcol.push_back(ntpur);
1175  nteffcol.push_back(nteff);
1176 
1177  ntnavmotherpdgcol.push_back(ntnavmotherpdg);
1178  ntnavmotherdangcol.push_back(ntnavmotherdang);
1179  ntnavmotheridcol.push_back(ntnavmotherid);
1180  ntnavmotherpcol.push_back(ntnavmotherp);
1181  ntnavmothervtxcol.push_back(ntnavmothervtx);
1182  ntnavmotherstopcol.push_back(ntnavmotherstop);
1183  ntnavmotherpurcol.push_back(ntnavmotherpur);
1184  ntnavmothereffcol.push_back(ntnavmothereff);
1185 
1186 
1187  }
1188 
1189 
1190  // return;//cell only
1191  //**********************************************************
1192  // Fill showers information
1193  //**********************************************************
1194 
1195  Vp4 showerpGam;
1196  showerpGam.clear();
1197 
1198  double showeretot = 0;
1199 
1200  TLorentzVector sh1p4(0,0,0,0);
1201  TLorentzVector trackpTot(0,0,0,0);
1202  TLorentzVector showerSumP4(0,0,0,0);
1203 
1204  TVector3 nuDir = geom->NuMIBeamDirection(); //(-0.845e-3, 52.552e-3, 0.998618);
1205  for(int i = 0; i < (int)slidcol.size(); i++){
1206  cet::maybe_ref<art::Ptr<rb::Shower> const> rshw(fos.at(i));
1207  art::Ptr<rb::Shower> shw = rshw.ref();
1208  // Set 4-Vector momentum
1209  double showerEnergy = fNuEEnergyAlg->ShowerEnergy(shw.get(), showercol, evt);
1210  double showereraw = showerEnergy;
1211  TVector3 p3(shw->Dir()[0]*showereraw,shw->Dir()[1]*showereraw,shw->Dir()[2]*showereraw);
1212  double maxll = -9999;
1213  int lltype = -1;
1214  std::map<int, float> tempLLL;
1215  std::map<int, float> tempLLT;
1216  tempLLL.clear();
1217  tempLLT.clear();
1218  tempLLL = slidcol[i]->fPartLongLL;
1219  tempLLT = slidcol[i]->fPartTransLL;
1220  for( int itype = 0; itype <11; ++itype){
1221  if(tempLLL[itype] + tempLLT[itype] > maxll){
1222  lltype = itype;
1223  maxll = tempLLL[itype] + tempLLT[itype];
1224  }
1225  }
1226  double m = xmass[0];
1227  if(lltype>=0&&lltype<7) m=xmass[lltype];
1228  double showerEtot0 = sqrt(p3.X()*p3.X()+p3.Y()*p3.Y()+p3.Z()*p3.Z());
1229  double P3 = sqrt(showerEtot0*showerEtot0+2*showerEtot0*m);
1230  double E = showerEtot0+m;
1231  TLorentzVector showerptrk(shw->Dir()[0]*P3,shw->Dir()[1]*P3,shw->Dir()[2]*P3,E);
1232  showerpGam.push_back(showerptrk);
1233  showeretot+= showerEnergy;
1234  showerSumP4+=showerptrk;
1235  }
1236  double evtSh1PlaneDedx[200];
1237  for(int plane=0;plane<200;plane++){
1238  evtSh1PlaneDedx[plane] = 0;
1239  }
1240 
1241  double evtSh1CellPlaneDedx[66];
1242  for(int cell=0;cell<11;cell++){
1243  for(int plane=0;plane<6;plane++){
1244  int cellplane = plane*11+cell;
1245  evtSh1CellPlaneDedx[cellplane] = 0;
1246  }
1247  }
1248 
1249  for(unsigned int i = 0; i < slidcol.size(); i++){
1250  cet::maybe_ref<art::Ptr<rb::Shower> const> rshw(fos.at(i));
1251  art::Ptr<rb::Shower> shw = rshw.ref();
1252 
1253  fParticleIDAlg->SetShower(shw.get(), showercol, evt);
1254 
1255  // General information
1256  double showerEnergy = fNuEEnergyAlg->ShowerEnergy(shw.get(), showercol, evt);
1257 
1258  int showerIsFidAna = 0;
1259  if ((fabs(shw->Start()[0]) < fXThresh) &&
1260  (fabs(shw->Stop()[0]) < fXThresh) &&
1261  (fabs(shw->Start()[1]) < fYThresh) &&
1262  (fabs(shw->Stop()[1]) < fYThresh) &&
1263  (std::min(shw->Start()[2],shw->Stop()[2]) > fZMin) &&
1264  (std::max(shw->Start()[2],shw->Stop()[2]) < fZMax)) showerIsFidAna = 1;
1265 
1266  int showerIsFid = 1;
1267  for (unsigned int pln=0; pln<shw->ExtentPlane(); ++pln){
1268  double plnrad = fParticleIDAlg->PlaneRadius(pln);
1269  TVector3 plncoord = fParticleIDAlg->PlaneHitXYZ(pln);
1270  double xdistedge = std::min(fabs(geom->DetHalfWidth()-plncoord[0]),fabs(-geom->DetHalfWidth()-plncoord[0]));
1271  double ydistedge = std::min(fabs(geom->DetHalfHeight()-plncoord[1]),fabs(-geom->DetHalfHeight()-plncoord[1]));
1272  double zdistedge = std::min(plncoord[2],fabs(geom->DetLength()-plncoord[2]));
1273  if ((xdistedge > 2.0*plnrad) && (ydistedge > 2.0*plnrad) && (zdistedge > 2.0*plnrad)) continue;
1274  showerIsFid = 0;
1275  break; //can stop checking once we have found one plane outside volume
1276  }
1277 
1278  m_showerTrueNuCCNC = evtTrueNuCCNC;
1279  m_showerTrueNuMode = evtTrueNuMode;
1280  m_showerTrueNuPdg = evtTrueNuPdg;
1281  m_showerTrueDang = ntdangcol[i];
1282  m_showerTruePdg = ntpdgcol[i];
1283  m_showerTrueP4[0] = ntpcol[i].Px();
1284  m_showerTrueP4[1] = ntpcol[i].Py();
1285  m_showerTrueP4[2] = ntpcol[i].Pz();
1286  m_showerTrueP4[3] = ntpcol[i].E();
1287 
1288  m_showerTrueNavMotherPdg = ntnavmotherpdgcol[i];
1289  m_showerTrueNavMotherDang = ntnavmotherdangcol[i];
1290  m_showerTrueNavMotherId = ntnavmotheridcol[i];
1291  m_showerTrueNavMotherP4[0] = ntnavmotherpcol[i].Px();
1292  m_showerTrueNavMotherP4[1] = ntnavmotherpcol[i].Py();
1293  m_showerTrueNavMotherP4[2] = ntnavmotherpcol[i].Pz();
1294  m_showerTrueNavMotherP4[3] = ntnavmotherpcol[i].E();
1295  m_showerTrueNavMotherVtx[0] = ntnavmothervtxcol[i][0];
1296  m_showerTrueNavMotherVtx[1] = ntnavmothervtxcol[i][1];
1297  m_showerTrueNavMotherVtx[2] = ntnavmothervtxcol[i][2];
1298  m_showerTrueNavMotherStop[0] = ntnavmotherstopcol[i][0];
1299  m_showerTrueNavMotherStop[1] = ntnavmotherstopcol[i][1];
1300  m_showerTrueNavMotherStop[2] = ntnavmotherstopcol[i][2];
1301 
1302 
1303  m_showerNPlane = shw->ExtentPlane();
1305  m_showerNTrk = slidcol.size();
1306  m_showerDedxANN = slidcol[i]->Value();
1307  m_showerDedxANNEPi0= slidcol[i]->ValueEPi0();
1308  m_showerDedxANNECos= slidcol[i]->ValueECos();
1309  if((int)i==ish1){
1310  m_showerIsSh1 = 1;
1311  }else if((int)i==ish2){
1312  m_showerIsSh1 = 2;
1313  }else{
1314  m_showerIsSh1 = 0;
1315  }
1316  m_showerStart[0] = (shw->Start())[0];
1317  m_showerStart[1] = (shw->Start())[1];
1318  m_showerStart[2] = (shw->Start())[2];
1319  m_showerStop[0] = (shw->Stop())[0];
1320  m_showerStop[1] = (shw->Stop())[1];
1321  m_showerStop[2] = (shw->Stop())[2];
1322  m_showerDir[0] = (shw->Dir())[0];
1323  m_showerDir[1] = (shw->Dir())[1];
1324  m_showerDir[2] = (shw->Dir())[2];
1325  m_showerIsFidAna = showerIsFidAna;
1327 
1328  // Fill longitudinal and transverse information
1329  for(int itd=0;itd<20;itd++){
1330  m_showerTransCellDedx[itd] = 0;
1331  }
1332 
1333  for(unsigned int plane=0;plane<200;plane++){
1334  m_showerPlaneDedx[plane] = 0;
1335  }
1336  for(unsigned int plane=0;plane<shw->ExtentPlane();plane++){
1337  if(plane>=200)break;
1338  double tdedx = fParticleIDAlg->PlaneLongDedx(plane);
1339  if (tdedx< 0) tdedx=0;
1340  m_showerPlaneDedx[plane] = tdedx;
1341  if((int)i==ish1) evtSh1PlaneDedx[plane] = tdedx;
1342  }
1343  for(int cell=0;cell<11;cell++){
1344  for(int plane=0;plane<6;plane++){
1345  int cellplane = plane*11+cell;
1346  double cpdedx = fParticleIDAlg->CellPlaneDedx(cell-5,plane);
1347  m_showerCellPlaneDedx[cellplane] = cpdedx;
1348  if((int)i==ish1) {
1349 // std::cout<<"cellplane, cpdedx "<<cellplane<<" "<<cpdedx<<std::endl;
1350  evtSh1CellPlaneDedx[cellplane] = cpdedx;
1351  }
1352  }
1353  }
1354 
1355  for(unsigned itp = 0;itp<20;itp++){
1356  double tdedx = fParticleIDAlg->PlaneTransDedx(itp);
1357  if (tdedx< 0) tdedx=0;
1358  m_showerTransCellDedx[itp] = tdedx;
1359  }
1360 
1361  m_showerRun= run;
1362  m_showerSubRun= srun;
1363  m_showerEvent= event;
1364  m_showerSubEvent= (int)sliceIdx;
1365  m_showerId= i;
1366  std::cout<<"GOING TO FILL SHOWER TREE\n";
1367  fShower->Fill();
1368  }
1369 
1370  //****************************************
1371  // Fill event level information
1372  //****************************************
1373 
1374 
1375  double mintopStart = 9999;
1376  double minbotStart = 9999;
1377  double mineastStart = 9999;
1378  double minwestStart = 9999;
1379  double minfrontStart = 9999;
1380  double minbackStart = 9999;
1381 
1382  double mintopStop = 9999;
1383  double minbotStop = 9999;
1384  double mineastStop = 9999;
1385  double minwestStop = 9999;
1386  double minfrontStop = 9999;
1387  double minbackStop = 9999;
1388  for(int i = 0; i < (int)slidcol.size(); i++){
1389  cet::maybe_ref<art::Ptr<rb::Shower> const> rshw(fos.at(i));
1390  art::Ptr<rb::Shower> shw = rshw.ref();
1391 
1392  TVector3 start = shw->Start();
1393  TVector3 stop = shw->Stop();
1394  mintopStart = std::min(mintopStart, livegeo->DistToTop(start));
1395  mintopStop = std::min(mintopStop, livegeo->DistToTop(stop ));
1396 
1397  minbotStart = std::min(minbotStart, livegeo->DistToBottom(start));
1398  minbotStop = std::min(minbotStop, livegeo->DistToBottom(stop ));
1399 
1400  mineastStart = std::min(mineastStart, livegeo->DistToEdgeXMinus(start));
1401  mineastStop = std::min(mineastStop, livegeo->DistToEdgeXMinus(stop ));
1402 
1403  minwestStart = std::min(minwestStart, livegeo->DistToEdgeXPlus(start));
1404  minwestStop = std::min(minwestStop, livegeo->DistToEdgeXPlus(stop ));
1405 
1406  minfrontStart = std::min(minfrontStart, livegeo->DistToFront(start));
1407  minfrontStop = std::min(minfrontStop, livegeo->DistToFront(stop ));
1408 
1409  minbackStart = std::min(minbackStart, livegeo->DistToBack(start));
1410  minbackStop = std::min(minbackStop, livegeo->DistToBack(stop ));
1411  }
1412  m_evtMinTop = std::min(mintopStart, mintopStop);
1413  m_evtMinBottom = std::min(minbotStart, minbotStop);
1414 
1415  m_evtMinFront = std::min(minfrontStart, minfrontStop);
1416  m_evtMinBack = std::min(minbackStart, minbackStop);
1417 
1418  m_evtMinEast = std::min(mineastStart, mineastStop);
1419  m_evtMinWest = std::min(minwestStart, minwestStop);
1420 
1421  if(ish1!=-1){
1422  cet::maybe_ref<art::Ptr<rb::Shower> const> rshw(fos.at(ish1));
1423  art::Ptr<rb::Shower> shw = rshw.ref();
1424 
1425  fParticleIDAlg->SetShower(shw.get(), showercol, evt);
1426 
1427  m_evtSh1Energy = fNuEEnergyAlg->ShowerEnergy(shw.get(), showercol, evt);
1428  if ((fabs(shw->Start()[0]) < fXThresh) &&
1429  (fabs(shw->Stop()[0]) < fXThresh) &&
1430  (fabs(shw->Start()[1]) < fYThresh) &&
1431  (fabs(shw->Stop()[1]) < fYThresh) &&
1432  (std::min(shw->Start()[2],shw->Stop()[2]) > fZMin) &&
1433  (std::max(shw->Start()[2],shw->Stop()[2]) < fZMax)) m_evtSh1IsFidAna = 1;
1434 
1435  m_evtSh1IsFid = 1;
1436  for (unsigned int pln=0; pln<shw->ExtentPlane(); ++pln){
1437  double plnrad = fParticleIDAlg->PlaneRadius(pln);
1438  TVector3 plncoord = fParticleIDAlg->PlaneHitXYZ(pln);
1439  double xdistedge = std::min(fabs(geom->DetHalfWidth()-plncoord[0]),fabs(-geom->DetHalfWidth()-plncoord[0]));
1440  double ydistedge = std::min(fabs(geom->DetHalfHeight()-plncoord[1]),fabs(-geom->DetHalfHeight()-plncoord[1]));
1441  double zdistedge = std::min(plncoord[2],fabs(geom->DetLength()-plncoord[2]));
1442  if ((xdistedge > 2.0*plnrad) && (ydistedge > 2.0*plnrad) && (zdistedge > 2.0*plnrad)) continue;
1443  m_evtSh1IsFid = 0;
1444  break; //can stop checking once we have found one plane outside the volume
1445  }
1446 
1447  m_evtSh1P4[0] = showerpGam[ish1].Px();
1448  m_evtSh1P4[1] = showerpGam[ish1].Py();
1449  m_evtSh1P4[2] = showerpGam[ish1].Pz();
1450  m_evtSh1P4[3] = showerpGam[ish1].E();
1451  m_evtSh1Start[0] = shw->Start()[0];
1452  m_evtSh1Start[1] = shw->Start()[1];
1453  m_evtSh1Start[2] = shw->Start()[2];
1454  m_evtSh1Stop[0] = shw->Stop()[0];
1455  m_evtSh1Stop[1] = shw->Stop()[1];
1456  m_evtSh1Stop[2] = shw->Stop()[2];
1457 
1458  m_evtSh1TrueDang = TMath::RadToDeg()*showerpGam[ish1].Angle(ntpcol[ish1].Vect());//ntdangcol[ish1];
1459  m_evtSh1TrueEDang = TMath::RadToDeg()*showerpGam[ish1].Angle(electruep4.Vect());
1460  m_evtSh1TruePdg = ntpdgcol[ish1];
1461  m_evtSh1TrueP4[0] = ntpcol[ish1].Px();
1462  m_evtSh1TrueP4[1] = ntpcol[ish1].Py();
1463  m_evtSh1TrueP4[2] = ntpcol[ish1].Pz();
1464  m_evtSh1TrueP4[3] = ntpcol[ish1].E();
1465  m_evtSh1TruePur = ntpurcol[ish1];
1466  m_evtSh1TrueEff = nteffcol[ish1];
1467  m_evtSh1TrueNavMotherPdg = ntnavmotherpdgcol[ish1];
1468  m_evtSh1TrueNavMotherDang = ntnavmotherdangcol[ish1];
1469  m_evtSh1TrueNavMotherId = ntnavmotheridcol[ish1];
1470  m_evtSh1TrueNavMotherP4[0] = ntnavmotherpcol[ish1].Px();
1471  m_evtSh1TrueNavMotherP4[1] = ntnavmotherpcol[ish1].Py();
1472  m_evtSh1TrueNavMotherP4[2] = ntnavmotherpcol[ish1].Pz();
1473  m_evtSh1TrueNavMotherP4[3] = ntnavmotherpcol[ish1].E();
1474  m_evtSh1TrueNavMotherPur = ntnavmotherpurcol[ish1];
1475  m_evtSh1TrueNavMotherEff = ntnavmothereffcol[ish1];
1476 
1477 
1478 
1479  double ellmax = -9999;
1480  double elllmax = -9999;
1481  double elltmax = -9999;
1482  std::map<int, float> tempLLL;
1483  std::map<int, float> tempLLT;
1484  tempLLL.clear();
1485  tempLLT.clear();
1486  tempLLL = slidcol[ish1]->fPartLongLL;
1487  tempLLT = slidcol[ish1]->fPartTransLL;
1488  for(int itype=0; itype != int(slid::DedxParticleType::kDEDXPARTICLETYPESIZE); ++itype){
1489  m_evtSh1DedxLLL[itype] = tempLLL[itype];
1490  m_evtSh1DedxLLT[itype] = tempLLT[itype];
1494  {
1495  if(ellmax<tempLLL[itype]+tempLLT[itype]){
1496  ellmax = tempLLL[itype]+tempLLT[itype];
1497  elllmax = tempLLL[itype];
1498  elltmax = tempLLT[itype];
1499  }
1500  }
1501  }
1502  m_evtSh1DedxLLL[0] = elllmax;
1503  m_evtSh1DedxLLT[0] = elltmax;
1504  for(int plane=0;plane<200;plane++){
1505  m_evtSh1PlaneDedx[plane] = evtSh1PlaneDedx[plane];
1506  }
1507  for(int cell=0;cell<11;cell++){
1508  for(int plane=0;plane<6;plane++){
1509  int cellplane = plane*11+cell;
1510  m_evtSh1CellPlaneDedx[cellplane] = evtSh1CellPlaneDedx[cellplane];
1511  }
1512  }
1513  m_evtSh1VtxGeV = slidcol[ish1]->VertexEnergy();
1514  m_evtSh1Pi0Mgg = slidcol[ish1]->Pi0mass();
1515  m_evtSh1Gap = slidcol[ish1]->Gap();
1516  m_evtSh1NPlane = shw->ExtentPlane();
1517  m_evtSh1NMIPPlane = slidcol[ish1]->NMIPPlanes();
1518  m_evtSh1NueEnergy = slidcol[ish1]->NueEnergy();
1519  double showerSumTheta = showerSumP4.Angle(nuDir);
1520  double showerSumPt = showerSumP4.Vect().Mag()*sin(showerSumTheta);
1521  m_evtSh1SumPt = showerSumPt;
1522  m_evtSh1SumP = showerSumP4.Vect().Mag();
1523  m_evtSh1CosTheta = slidcol[ish1]->CosTheta();
1524  m_evtSh1Length = shw->TotalLength();
1525  m_evtSh1DedxANN = slidcol[ish1]->Value();
1526  m_evtSh1DedxANNEPi0 = slidcol[ish1]->ValueEPi0();
1527  m_evtSh1DedxANNECos = slidcol[ish1]->ValueECos();
1528 
1529  if(ish2!=-1){
1530  cet::maybe_ref<art::Ptr<rb::Shower> const> rshw2(fos.at(ish2));
1531  art::Ptr<rb::Shower> shw2 = rshw2.ref();
1532  m_evtSh2Energy = fNuEEnergyAlg->ShowerEnergy(shw.get(), showercol, evt);
1533 
1534  cet::maybe_ref<art::Ptr<rb::Shower> const> rshw1(fos.at(ish1));
1535  art::Ptr<rb::Shower> shw1 = rshw1.ref();
1536 
1537  m_evtPi0StartDist = (shw1->Start()-shw2->Start()).Mag();
1538 
1539  double x1 = shw1->Start()[0];
1540  double y1 = shw1->Start()[1];
1541  double z1 = shw1->Start()[2];
1542  double p1 = (shw1->Stop()[0] - shw1->Start()[0])/(shw1->Stop() - shw1->Start()).Mag();
1543  double q1 = (shw1->Stop()[1] - shw1->Start()[1])/(shw1->Stop() - shw1->Start()).Mag();
1544  double r1 = (shw1->Stop()[2] - shw1->Start()[2])/(shw1->Stop() - shw1->Start()).Mag();
1545  double x2 = shw2->Start()[0];
1546  double y2 = shw2->Start()[1];
1547  double z2 = shw2->Start()[2];
1548  double p2 = (shw2->Stop()[0] - shw2->Start()[0])/(shw2->Stop() - shw2->Start()).Mag();
1549  double q2 = (shw2->Stop()[1] - shw2->Start()[1])/(shw2->Stop() - shw2->Start()).Mag();
1550  double r2 = (shw2->Stop()[2] - shw2->Start()[2])/(shw2->Stop() - shw2->Start()).Mag();
1551 
1552  TMatrixD dNorm(3,3);
1553  dNorm[0][0] = x2-x1;
1554  dNorm[0][1] = y2-y1;
1555  dNorm[0][2] = z2-z1;
1556  dNorm[1][0] = p1;
1557  dNorm[1][1] = q1;
1558  dNorm[1][2] = r1;
1559  dNorm[2][0] = p2;
1560  dNorm[2][1] = q2;
1561  dNorm[2][2] = r2;
1562 
1563  double ma = x2-x1;
1564  double mb = y2-y1;
1565  double mc = z2-z1;
1566  double md = p1;
1567  double me = q1;
1568  double mf = r1;
1569  double mg = p2;
1570  double mh = q2;
1571  double mi = r2;
1572 
1573 
1574  TMatrixD dDeNorm1(2,2);
1575  TMatrixD dDeNorm2(2,2);
1576  TMatrixD dDeNorm3(2,2);
1577  dDeNorm1[0][0] = p1;
1578  dDeNorm1[0][1] = q1;
1579  dDeNorm1[1][0] = p2;
1580  dDeNorm1[1][1] = q2;
1581  double det1 = dDeNorm1[0][0]*dDeNorm1[1][1] - dDeNorm1[0][1]*dDeNorm1[1][0];
1582  dDeNorm2[0][0] = q1;
1583  dDeNorm2[0][1] = r1;
1584  dDeNorm2[1][0] = q2;
1585  dDeNorm2[1][1] = r2;
1586  double det2 = dDeNorm2[0][0]*dDeNorm2[1][1] - dDeNorm2[0][1]*dDeNorm2[1][0];
1587  dDeNorm3[0][0] = r1;
1588  dDeNorm3[0][1] = p1;
1589  dDeNorm3[1][0] = r2;
1590  dDeNorm3[1][1] = p2;
1591  double det3 = dDeNorm3[0][0]*dDeNorm3[1][1] - dDeNorm3[0][1]*dDeNorm3[1][0];
1592  double det0 = ma*me*mi + mb*mf*mg + mc*md*mh - mc*me*mg - mb*md*mi - ma*mf*mh;
1593  double pi0LineDist = fabs(det0)/sqrt(det1*det1+det2*det2+det3*det3);
1594  m_evtPi0LineDist = pi0LineDist;
1595  }
1596  }
1597  m_evtEtot = showeretot;
1598  fEvent->Fill();
1599  }
1600  return;
1601  }
1602 
1604  {
1605  std::cout<<"Total number of events "<<Ncount0<<std::endl;
1606  }
1607 
1609 
1610 }
1611 //......................................................................
NuEEnergyAlg * fNuEEnergyAlg
T max(const caf::Proxy< T > &a, T b)
virtual void reconfigure(const fhicl::ParameterSet &pset)
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
T1 fmod(T1 number1, T2 number2)
Definition: d0nt_math.hpp:102
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Cluster.cxx:121
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
Double_t angle
Definition: plot.C:86
virtual void beginSubRun(const art::SubRun &sr)
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
double m_showerTrueNavMotherVtx[3]
Trajectory class.
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
double EndZ() const
Definition: MCParticle.h:227
Float_t y1[n_points_granero]
Definition: compare.C:5
double m_evtSh1TruePiPlaneEnergy[10]
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
double m_showerCellPlaneDedx[66]
virtual void beginRun(const art::Run &run)
X or Y views.
Definition: PlaneGeo.h:30
virtual void analyze(const art::Event &evt)
double m_evtSh1TrueMuPlaneEnergy[10]
Float_t x1[n_points_granero]
Definition: compare.C:5
TVector3 PlaneHitXYZ(unsigned int pIdx)
Return point of intersection between shower core and plane.
std::vector< int > SliceToOrderedNuIdsByEff(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable) const
double m_showerTrueMuPlaneEnergy[10]
ParticleIDAlg * fParticleIDAlg
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
T sqrt(T number)
Definition: d0nt_math.hpp:156
double m_showerTrueNPlaneEnergy[10]
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
A collection of associated CellHits.
Definition: Cluster.h:47
double DetLength() const
std::vector< NeutrinoWithIndex > allNeutrinoInteractions() const
Function primarily for CAFMaker - returns a vector of all MCTruth interactions along with their truth...
double PlaneRadius(unsigned int pIdx)
return shower radius for a plane index
void abs(TH1 *hist)
static const int NUMUID
double m_evtSh1TruePPlaneEnergy[10]
A single unit of energy deposition in the liquid scintillator.
Definition: FLSHit.h:19
double DistToFront(TVector3 vertex)
std::vector< int > Vint
double DistToBack(TVector3 vertex)
bool isRealData() const
Definition: Event.h:83
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
double EndY() const
Definition: MCParticle.h:226
double PlaneLongDedx(unsigned int pIdx)
return longitudinal dedx for a specified plane index
static const int PIMINUSID
double m_showerTruePiPlaneEnergy[10]
object containing MC flux information
Definition: Run.h:31
double HitCollectionEfficiency(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, const std::vector< const rb::CellHit * > &allhits, const geo::View_t &view, std::map< int, double > *effMap=0, bool energyEff=false, double *desiredEnergy=0, double *totalEnergy=0, int *desiredHits=0, int *totalHits=0) const
Returns the fraction of all energy in an event from a specific set of Geant4 track IDs that are repre...
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
fhicl::ParameterSet fNuEEnergyAlgPSet
FCL parameters for NuE Energy alorithm.
unsigned int PlaneToGlobalPlane(unsigned int pIdx)
Return the gloabel plane index for an index in the shower coordinate.
virtual void beginJob()
int GetPDG() const
PDG.
Definition: FLSHit.h:43
Double_t q2[12][num]
Definition: f2_nu.C:137
double m_showerTruePlaneEnergy[10]
double DistToTop(TVector3 vertex)
TString part[npart]
Definition: Style.C:32
Definition: NueSkimmer.h:24
double DistToBottom(TVector3 vertex)
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
double CellPlaneDedx(int tpIdx, unsigned int pIdx)
return transverse dedx for specified transverse plane
double m_evtSh1PlaneDedx[200]
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
Float_t E
Definition: plot.C:20
double m_evtSh1TrueNPlaneEnergy[10]
const double xmass[7]
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
std::vector< sim::FLSHit > HitToFLSHit(const rb::CellHit &hit) const
All the FLSHits that contributed to this hit, sorted from most to least light.
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
float PE() const
Definition: CellHit.h:42
double m_showerTruePPlaneEnergy[10]
Collect Geo headers and supply basic geometry functions.
double m_evtSh1CellPlaneDedx[66]
caf::StandardRecord * sr
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
double m_showerPlaneDedx[200]
double m_showerTrueEPlaneEnergy[10]
double showerEnergy
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Vertex location in position and time.
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:308
double HitCollectionPurity(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, std::map< int, double > *purMap=0, std::map< int, int > *parents=0, bool energyPur=false) const
Returns the fraction of hits in a collection that come from the specified Geant4 track ids...
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
OStream cout
Definition: OStream.cxx:6
double ShowerEnergy(const rb::Shower *shower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Returns the total energy of a shower. along with corrections due to dead material and threshold effec...
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
static const int NUEID
double m_showerTrueGPlaneEnergy[10]
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
EventNumber_t event() const
Definition: EventID.h:116
std::vector< TrackIDE > HitsToTrackIDE(const std::vector< const rb::CellHit * > &hits, bool useBirksE=false) const
Returns vector of TrackIDE structs contributing to the given collection of hits.
fhicl::ParameterSet fParticleIDAlgPSet
FCL parameters for particle ID alorithm (loglikelihoods and dE/dx)
double m_evtSh1TruePlaneEnergy[10]
T * make(ARGS...args) const
std::vector< double > Vdouble
float GeV() const
Definition: RecoHit.cxx:69
T sin(T number)
Definition: d0nt_math.hpp:132
double PlaneTransDedx(unsigned int tpIdx)
return transverse dedx for specified transverse plane
Calculate dE/dx and log likelihoods for different particle hypotheses. This information will be of us...
double DetHalfWidth() const
double m_evtSh1TrueEPlaneEnergy[10]
double DistToEdgeXMinus(TVector3 vertex)
T const * get() const
Definition: Ptr.h:321
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
LIDTraining(fhicl::ParameterSet const &pset)
Calculates deposited and corrected energy of the electron shower and of electron flavoured neutrino...
Definition: structs.h:12
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
static const int NUTAUID
double m_showerTransCellDedx[20]
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
void geom(int which=0)
Definition: geom.C:163
float Mag() const
std::vector< TLorentzVector > Vp4
Build slid::LID objects to store electron ID, if asked for, otherwise, calculate LID info and make av...
Definition: FillPIDs.h:13
double showerIsFid
bool NeutrinoSet() const
Definition: MCTruth.h:77
float GetEdep() const
Get total Energy deposited into the cell for the whole FLSHit.
Definition: FLSHit.h:31
T min(const caf::Proxy< T > &a, T b)
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
double m_showerTrueNavMotherP4[4]
Int_t trackID
Definition: plot.C:84
double EndX() const
Definition: MCParticle.h:225
RunNumber_t run() const
Definition: Event.h:77
static constexpr Double_t mg
Definition: Munits.h:210
Event generator information.
Definition: MCNeutrino.h:18
int Mode() const
Definition: MCNeutrino.h:149
double m_showerTrueNavMotherStop[3]
std::vector< NeutrinoEffPur > SliceToMCTruth(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of MCTruth, efficiency, and purity of that neutrino interaction relative to the slice. Efficiency is defined as FLS energy from neutrino interaction in slice / total FLS energy from neutrino interaction in event. This vector is sorted from the highest efficiency interaction to lowest. This function returns all MCTruth, including those without neutrino interactions, i.e. cosmics.
double m_evtSh1TrueNavMotherP4[4]
double m_evtSh1TrueGPlaneEnergy[10]
EventID id() const
Definition: Event.h:56
static const int GAMMAID
Encapsulate the geometry of one entire detector (near, far, ndos)
std::string fPotLabel
Module that produced the POTSum object.
TVector3 Vect() const
double DistToEdgeXPlus(TVector3 vertex)
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const
void SetShower(const rb::Shower *vShower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Set rb::Prong to be analyzed. This must be set before any calculations can be done.
static const int PIPLUSID
int EveId(const int trackID) const
enum BeamMode string
static const int PIZEROID