RecoJMShower_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 /// \brief A producer that reconstruct shower energy, topological variables
3 /// and PID information
4 /// \author Jianming Bian - bianjm@umn.physics.edu
5 ////////////////////////////////////////////////////////////////////////
6 extern "C" {
7 #include <sys/types.h>
8 #include <sys/stat.h>
9 }
10 
11 #include <fstream>
12 #include <cstdlib>
13 #include <vector>
14 #include <string>
15 
16 // ROOT includes
17 #include "TFile.h"
18 #include "TH1D.h"
19 #include "TVector3.h"
20 #include "TTree.h"
21 #include "TMultiLayerPerceptron.h"
22 
23 // ART includes
24 #include "art/Framework/Core/EDProducer.h"
25 #include "art/Framework/Core/ModuleMacros.h"
26 #include "art/Framework/Principal/Event.h"
27 #include "art/Framework/Principal/Handle.h"
28 #include "art/Framework/Services/Registry/ServiceHandle.h"
29 #include "art/Framework/Services/Optional/TFileService.h"
30 //#include "art/Framework/Core/FindManyP.h"
31 //#include "art/Persistency/Common/Ptr.h"
32 #include "canvas/Persistency/Common/FindManyP.h"
33 #include "canvas/Persistency/Common/Ptr.h"
34 #include "cetlib/search_path.h"
35 
36 // NOvA includes
37 #include "Calibrator/Calibrator.h"
38 #include "CalibrationUtils/CalibUtil.h"
39 #include "Geometry/Geometry.h"
40 #include "GeometryObjects/CellGeo.h"
41 //#include "Geometry/func/CellUniqueId.h"
42 #include "RecoBase/CellHit.h"
43 #include "RecoBase/Cluster.h"
44 #include "RecoBase/FilterList.h"
45 #include "RecoBase/RecoHit.h"
46 #include "RecoBase/Track.h"
47 #include "RecoBase/Shower.h"
48 #include "RecoBase/Vertex.h"
49 #include "nusimdata/SimulationBase/MCParticle.h"
50 #include "nusimdata/SimulationBase/MCTruth.h"
51 #include "Simulation/FLSHitList.h"
52 #include "Utilities/func/EnvExpand.h"
53 #include "Utilities/func/MathUtil.h"
54 #include "Utilities/AssociationUtil.h"
55 #include "DAQChannelMap/DAQChannelMapConstants.h"
56 //#include "Geometry/OfflineChan.h"
57 #include "RecoJMShower/JMShower.h"
58 
59 
60 namespace jmshower {
61 
62  /// Store some RVP-like variables that could potentially factor into the EID
63  struct RVPInfo
64  {
65  double efrac_2slides; ///< highest energy fraction in 2 plane sliding window out of slice
66  double efrac_4slides; ///< highest energy fraction in 4 plane sliding window out of slice
67  double efrac_6slides; ///< highest energy fraction in 6 plane sliding window out of slice
68  double efrac_20slides; ///< highest energy fraction in first 20 planes
69  double efrac_2sigmaHOR; ///< energy fraction within 2sigma of slice center
70  double efrac_4sigmaHOR; ///< energy fraction within 4sigma of slice center
71 
72  double eiso_2sigmaHOR; ///< energy isolation of 2sigma raod
73  double eiso_4sigmaHOR; ///< energy isolation of 4sigma raod
74  double eiso_6sigmaHOR; ///< energy isolation of 6sigma raod
75  };
76 
77  class RecoJMShower : public art::EDProducer
78  {
79  public:
80  RecoJMShower(fhicl::ParameterSet const& pset);
81  ~RecoJMShower();
82 
83  void preEvent(art::Event const& evt);
84  void postBeginRun(art::Run const& run);
85  void reconfigure(const fhicl::ParameterSet &pset);
86  void produce(art::Event& evt);
87  void beginRun(art::Run& run);
88 
89  // Histograms for probability calculation
90  TH1D *fHtPlaneDedxE[4][11][200];
91  TH1D *fHtPlaneDedxG[4][11][200];
92  TH1D *fHtPlaneDedxMu[4][11][200];
93  TH1D *fHtPlaneDedxPi0[4][11][200];
94  TH1D *fHtPlaneDedxHad[4][11][200];
95  TH1D *fHtPlaneDedxP[4][11][200];
96  TH1D *fHtPlaneDedxN[4][11][200];
97  TH1D *fHtPlaneDedxPi[4][11][200];
98  TH1D *fHtPlaneDedxEqe[4][11][200];
99  TH1D *fHtPlaneDedxEres[4][11][200];
100  TH1D *fHtPlaneDedxEdis[4][11][200];
101  TH1D *fHtPlaneDedxEcoh[4][11][200];
102  TH1D *fHtPlaneDedxEsg[4][41][200];
103 
104  unsigned fHtExpPlaneE[4][11];
105  unsigned fHtExpPlaneG[4][11];
106  unsigned fHtExpPlaneMu[4][11];
107  unsigned fHtExpPlanePi0[4][11];
108  unsigned fHtExpPlaneHad[4][11];
109  unsigned fHtExpPlaneP[4][11];
110  unsigned fHtExpPlaneN[4][11];
111  unsigned fHtExpPlanePi[4][11];
112  unsigned fHtExpPlaneEqe[4][11];
113  unsigned fHtExpPlaneEres[4][11];
114  unsigned fHtExpPlaneEdis[4][11];
115  unsigned fHtExpPlaneEcoh[4][11];
116  unsigned fHtExpPlaneEsg[4][41];
117 
118 
119 
120 /*
121  TH1D *fHtmipPlaneDedxE[11][200];
122  TH1D *fHtmipPlaneDedxG[11][200];
123  TH1D *fHtmipPlaneDedxMu[11][200];
124  TH1D *fHtmipPlaneDedxPi0[11][200];
125  TH1D *fHtmipPlaneDedxHad[11][200];
126  TH1D *fHtmipPlaneDedxP[11][200];
127  TH1D *fHtmipPlaneDedxN[11][200];
128  TH1D *fHtmipPlaneDedxPi[11][200];
129  TH1D *fHtmipPlaneDedxEqe[11][200];
130  TH1D *fHtmipPlaneDedxEres[11][200];
131  TH1D *fHtmipPlaneDedxEdis[11][200];
132  TH1D *fHtmipPlaneDedxEcoh[11][200];
133 
134  TH1D *fHtshowerPlaneDedxE[11][200];
135  TH1D *fHtshowerPlaneDedxG[11][200];
136  TH1D *fHtshowerPlaneDedxMu[11][200];
137  TH1D *fHtshowerPlaneDedxPi0[11][200];
138  TH1D *fHtshowerPlaneDedxHad[11][200];
139  TH1D *fHtshowerPlaneDedxP[11][200];
140  TH1D *fHtshowerPlaneDedxN[11][200];
141  TH1D *fHtshowerPlaneDedxPi[11][200];
142  TH1D *fHtshowerPlaneDedxEqe[11][200];
143  TH1D *fHtshowerPlaneDedxEres[11][200];
144  TH1D *fHtshowerPlaneDedxEdis[11][200];
145  TH1D *fHtshowerPlaneDedxEcoh[11][200];
146 
147  TH2F *fHt2dPlaneDedxE[11][200];
148  TH2F *fHt2dPlaneDedxG[11][200];
149  TH2F *fHt2dPlaneDedxMu[11][200];
150  TH2F *fHt2dPlaneDedxPi0[11][200];
151  TH2F *fHt2dPlaneDedxHad[11][200];
152  TH2F *fHt2dPlaneDedxP[11][200];
153  TH2F *fHt2dPlaneDedxN[11][200];
154  TH2F *fHt2dPlaneDedxPi[11][200];
155  TH2F *fHt2dPlaneDedxEqe[11][200];
156  TH2F *fHt2dPlaneDedxEres[11][200];
157  TH2F *fHt2dPlaneDedxEdis[11][200];
158  TH2F *fHt2dPlaneDedxEcoh[11][200];
159 */
160 
161  TH1D *fHtCellTransDedxE[4][11][20];
162  TH1D *fHtCellTransDedxG[4][11][20];
163  TH1D *fHtCellTransDedxMu[4][11][20];
164  TH1D *fHtCellTransDedxPi0[4][11][20];
165  TH1D *fHtCellTransDedxHad[4][11][20];
166  TH1D *fHtCellTransDedxP[4][11][20];
167  TH1D *fHtCellTransDedxN[4][11][20];
168  TH1D *fHtCellTransDedxPi[4][11][20];
169  TH1D *fHtCellTransDedxEqe[4][11][20];
170  TH1D *fHtCellTransDedxEres[4][11][20];
171  TH1D *fHtCellTransDedxEdis[4][11][20];
172  TH1D *fHtCellTransDedxEcoh[4][11][20];
173  TH1D *fHtCellTransDedxEsg[4][41][20];
174 
175 
176 
177  std::string fFnameWeightsShape;
178  std::string fFnameWeightsShapeElec;
179  std::string fFnameWeightsShapeQE;
180  std::string fFnameWeightsShapeRES;
181  std::string fFnameWeightsShapeDIS;
182  TMultiLayerPerceptron *mlp_shape;
183  TMultiLayerPerceptron *mlp_shape_elec;
184  TMultiLayerPerceptron *mlp_shape_E;
185  TMultiLayerPerceptron *mlp_shape_NoCos;
186  TMultiLayerPerceptron *mlp_shape_E_NoCos;
187  TMultiLayerPerceptron *mlp_shape_QE;
188  TMultiLayerPerceptron *mlp_shape_RES;
189  TMultiLayerPerceptron *mlp_shape_DIS;
190  TMultiLayerPerceptron *mlp_shape_EL;
191  TMultiLayerPerceptron *mlp_shape_epi0;
192  TMultiLayerPerceptron *mlp_shape_epi0EL;
193 
194 // Energy deconvolution map
195  std::map<geo::OfflineChan, std::vector<double> > fHitEWeightMap;
196  std::map<geo::OfflineChan, std::vector<double> > fHitDistMap;
197  std::map<geo::OfflineChan, double> fHitEMap;
198  std::vector<std::vector<jmshower::PCluster> > fShowerPClusterCol;
199  std::vector<std::vector<jmshower::PCluster> > fTrackPClusterCol;
200  std::vector<std::vector<jmshower::TCCluster> > fShowerTCClusterCol;
201  std::vector<jmshower::JMShower> fExcludeShowerCol;
202 
203 
204 
205 // Geometry
206  double GetPointDoca(TVector3 vtx, TVector3 start, TVector3 stop);
207  TVector3 GetCellDistToPoint(unsigned int plane, unsigned int cell, TVector3 point);
208  double GetCellDistToTrk(unsigned int plane, unsigned int cell, TVector3 start, TVector3 stop);
209 
210  TVector3 GetShwStop(rb::Prong shw);//Calculate the end point of a shower
211  TVector3 GetShwStop(art::Ptr<rb::Prong> shw);//Calculate the end point of a shower
212 
213 
214  TVector3 GetTrkHitPos(unsigned int plane, unsigned int cell, art::Ptr<rb::Prong> trk);
215  TVector3 GetTrkHitPos(unsigned int plane, unsigned int cell, rb::Prong trk);
216  TVector3 GetTrkHitPos(unsigned int plane, unsigned int cell, TVector3 trkdir, TVector3 trkstart, TVector3 trkstop);
217  TVector3 GetTrkPlanePos(unsigned int plane, rb::Prong trk);
218  unsigned int GetTrkPlaneCell(unsigned int plane, art::Ptr<rb::Prong> trk);
219  unsigned int GetTrkPlaneCell(unsigned int plane, rb::Prong trk);
220  unsigned int GetTrkCPlaneCell(unsigned int plane, unsigned int i);
221  TVector3 GetTrkPlaneDistToEdge(unsigned int plane, art::Ptr<rb::Prong> trk);
222  TVector3 GetTrkPlaneDistToEdge(unsigned int plane, rb::Prong trk);
223  bool IsFiducial(rb::Prong shower);
224  double GetTrkDistToEdge(rb::Prong shower);
225  double GetTrkHitPath(unsigned int plane, unsigned int cell, rb::Prong shower);
226  unsigned int ContStartPlane(rb::Prong shower, unsigned int startp, int contp);
227 // Energy
228  double DepositEnergy(rb::Cluster cluster, int i);
229  double DepositEnergy(rb::Prong shower);
230  double Energy(rb::Cluster cluster, int i);
231  double Energy(rb::Prong shower);
232  double ECorrMC(double gev);
233  double CellADCToGeV(double d, geo::View_t view, int detid);
234  double CellEDataMC(double gev, double w, geo::View_t view, int detid); // correct cell DATA/MC difference
235 
236  double PlaneDedxDataMC(double gev, int detid); // correct plane dedx DATA/MC difference
237  double CellTransDataMC(double gev, int detid); // correct transverse dedx DATA/MC difference
238 
239  double Radius(rb::Prong shower);
240  double Radius(rb::Cluster cluster, rb::Prong shower);
241  double RadiusV(rb::Cluster cluster);
242 
243 // Shower shape
244  double Gap(rb::Prong shower, TVector3 vertex);
245  double GetPlaneDedx(rb::Prong shower, unsigned int plane);
246  double GetPlaneDedxProb(rb::Prong shower, int ireg, int igev, unsigned int plane, double dedx, const int type);
247 // double GetPlaneDedxProb(rb::Prong shower, int igev, unsigned int plane, double dedx, unsigned int nmip, const int type);
248  double GetInterPlaneDedxProb(rb::Prong shower, unsigned int plane, double dedx, const int type, TVector3 pos, int iReg, int igev0, int igev1, double e0, double e1);
249 // double GetInterPlaneDedxProb(rb::Prong shower, unsigned int plane, double dedx, const int type, unsigned int nmip);
250  double GetDedxLongLL(rb::Prong shower, const int type);
251  double GetDedxInvLongLL(rb::Prong shower, const int type);
252  double GetDedxLongLL(rb::Prong shower, const int type, unsigned int startp, int contp);
253 
254  double GetCellTransDedx(rb::Prong shower, unsigned int cell);
255  double GetCellTransDedxProb(rb::Prong shower, int ireg, int igev, unsigned int cell, double dedx, const int type);
256  double GetInterCellTransDedxProb(rb::Prong shower, unsigned int cell, double dedx, const int type, TVector3 pos, int iReg, int igev0, int igev1, double e0, double e1);
257  double GetDedxTransLL(rb::Prong shower, const int type);
258  bool IsMuon(jmshower::JMShower shower);
259  bool IsInvPhoton(jmshower::JMShower shower);
260 
261  unsigned int NMIPPlane(rb::Prong shower);
262  unsigned int NMIPPlane(rb::Prong shower, unsigned int startp);
263  double CalcANN(unsigned int mode, jmshower::JMShower shower);
264 
265  //function to store rvp type variables that could potentially be used
266  RVPInfo GetRVPStats(art::Handle<std::vector<rb::Cluster> > slices, unsigned int ic);
267  TVector3 GetCellNodePos(unsigned int plane1, unsigned int cell1, unsigned int plane2, unsigned int cell2);
268 
269  std::string fLibPath;
270  std::string pSliceDir;
271  std::string pVertexLabel;
272  std::string pInstLabel;
273  std::string pTrackDir;
274  double pCellHitPHThresh;
275  int pPIDHist;
276  int pRecoANN;
277  bool pRecoPID;
278  bool pRecoRVP;
279  bool pRecoNode;
280  bool pRecoPlaneR;
281  bool pRecoInvLL;
282  bool pRecluster;
283  bool pRecoVtx;
284  bool pVtxFit;
285  bool pBadChannel;
286  bool pCorAbsCell;
287  bool pCorPigtail;
288  bool pCorDataCell;
289  bool pCorDataDedx;
290  bool pUseUncalib;
291 
292  int pXEdgeId1;
293  int pXEdgeId2;
294  int pXEdgeId3;
295  int pXEdgeId4;
296  int pYEdgeId1;
297  int pYEdgeId2;
298  int pYEdgeId3;
299  int pYEdgeId4;
300  int pZEdgeId1;
301  int pZEdgeId2;
302  int pZEdgeId3;
303  int pZEdgeId4;
304 
305  double pXEdge1;
306  double pXEdge2;
307  double pXEdge3;
308  double pXEdge4;
309  double pYEdge1;
310  double pYEdge2;
311  double pYEdge3;
312  double pYEdge4;
313  double pZEdge1;
314  double pZEdge2;
315  double pZEdge3;
316  double pZEdge4;
317 
318 // Exclude cluster
319  TVector3 GetCentroid(rb::Cluster cluster);
320  int GetPlaneE1Cell(jmshower::PCluster pcluster);
321  int GetPlaneCentroidCell(jmshower::PCluster pcluster, rb::Prong shower);
322 
323  public:
324  std::vector<jmshower::JMShower> RecoShowers(art::EDProducer const& prod,
325  art::Event& evt,
326  std::unique_ptr< std::vector<jmshower::JMShower> >& recoshowercol,
327  std::unique_ptr< art::Assns<jmshower::JMShower, rb::Prong> >&);
328 
329 // std::vector<jmshower::JMShower> RecoShowers(const art::Event& evt);
330  std::vector<jmshower::JMShower> ExclShowers(const art::Event& evt);
331 
332  private:
333  bool PreSel(const art::Event& evt);
334  bool LoadDedxHistograms(int detid);
335  bool LoadTreeWeights(int detid);
336 
337  int fDetid;
338  bool isrealdata;
339  bool fLoaded; ///< Have we done beginRun once yet?
340  std::vector<int> fprongidx;
341  };
342 
343 } // end namespace jmshower
344 
345 namespace jmshower
346 {
347 //------------------------------------------------------------
348  PCluster::PCluster() : fP(9999)
349  {
350  }
351 
352  PCluster::~PCluster()
353  {
354  }
355 
356  TCCluster::TCCluster() : fTc(9999)
357  {
358  }
359 
360  TCCluster::~TCCluster()
361  {
362  }
363 
364 
365 //------------------------------------------------------------
366  RecoJMShower::RecoJMShower(fhicl::ParameterSet const& pset):
367  pSliceDir(pset.get<std::string>("SliceDir")),
368  pTrackDir(pset.get<std::string>("TrackDir")),
369  pCellHitPHThresh(50.),
370  pPIDHist(1),
371  pRecoANN(0),
372  pRecoPID(true),
373  pRecoRVP(true),
374  pRecoNode(true),
375  pRecoPlaneR(true),
376  pRecoInvLL(true),
377  pRecluster(true),
378  pRecoVtx(true),
379  pVtxFit(false),
380  pBadChannel(false),
381  pCorAbsCell(false),
382  pCorPigtail(true),
383  pCorDataCell(false),
384  pCorDataDedx(false),
385  pUseUncalib(false),
386  pXEdge1(-774.7),
387  pXEdge2(774.7),
388  pXEdge3(-774.7),
389  pXEdge4(-774.7),
390  pYEdge1(-774.7),
391  pYEdge2(774.7),
392  pYEdge3(-774.7),
393  pYEdge4(-774.7),
394  pZEdge1(0.0),
395  pZEdge2(5959.08),
396  pZEdge3(0.0),
397  pZEdge4(0.0),
398  fLoaded(false)
399  {
400  reconfigure(pset);
401  produces< std::vector<jmshower::JMShower> >();
402  produces< art::Assns<jmshower::JMShower, rb::Cluster> >();
403  produces< art::Assns<jmshower::JMShower, rb::Prong> >();
404  }
405 
406  //------------------------------------------------------------
407  RecoJMShower::~RecoJMShower()
408  {
409  }
410 
411  //------------------------------------------------------------
412  void RecoJMShower::reconfigure(const fhicl::ParameterSet& pset)
413  {
414  fLibPath = pset.get< std::string >("LibraryPath");
415  pSliceDir = pset.get< std::string >("SliceDir");
416  pVertexLabel = pset.get< std::string >("VertexLabel");
417  pTrackDir = pset.get< std::string >("TrackDir");
418  pInstLabel = pset.get< std::string >("InstLabel");
419  pCellHitPHThresh = pset.get< double >("CellHitPHThresh");
420  pPIDHist = pset.get< int >("PIDHist");
421  pRecoANN = pset.get< int >("RecoANN");
422  pRecoPID = pset.get< bool >("RecoPID");
423  pRecoRVP = pset.get< bool >("RecoRVP");
424  pRecoNode = pset.get< bool >("RecoNode");
425  pRecoPlaneR = pset.get< bool >("RecoPlaneR");
426  pRecoInvLL = pset.get< bool >("RecoInvLL");
427  pRecluster = pset.get< bool >("Recluster");
428  pRecoVtx = pset.get< bool >("RecoVtx");
429  pVtxFit = pset.get< bool >("VtxFit");
430  pBadChannel = pset.get< bool >("BadChannel");
431  pCorAbsCell = pset.get< bool >("CorAbsCell");
432  pCorPigtail = pset.get< bool >("CorPigtail");
433  pCorDataCell = pset.get< bool >("CorDataCell");
434  pCorDataDedx = pset.get< bool >("CorDataDedx");
435  pUseUncalib = pset.get< bool >("UseUncalib");
436  pXEdgeId1 = pset.get< int >("XEdgeId1");
437  pXEdgeId2 = pset.get< int >("XEdgeId2");
438  pXEdgeId3 = pset.get< int >("XEdgeId3");
439  pXEdgeId4 = pset.get< int >("XEdgeId4");
440  pYEdgeId1 = pset.get< int >("YEdgeId1");
441  pYEdgeId2 = pset.get< int >("YEdgeId2");
442  pYEdgeId3 = pset.get< int >("YEdgeId3");
443  pYEdgeId4 = pset.get< int >("YEdgeId4");
444  pZEdgeId1 = pset.get< int >("ZEdgeId1");
445  pZEdgeId2 = pset.get< int >("ZEdgeId2");
446  pZEdgeId3 = pset.get< int >("ZEdgeId3");
447  pZEdgeId4 = pset.get< int >("ZEdgeId4");
448  pXEdge1 = pset.get< double >("XEdge1");
449  pXEdge2 = pset.get< double >("XEdge2");
450  pXEdge3 = pset.get< double >("XEdge3");
451  pXEdge4 = pset.get< double >("XEdge4");
452  pYEdge1 = pset.get< double >("YEdge1");
453  pYEdge2 = pset.get< double >("YEdge2");
454  pYEdge3 = pset.get< double >("YEdge3");
455  pYEdge4 = pset.get< double >("YEdge4");
456  pZEdge1 = pset.get< double >("ZEdge1");
457  pZEdge2 = pset.get< double >("ZEdge2");
458  pZEdge3 = pset.get< double >("ZEdge3");
459  pZEdge4 = pset.get< double >("ZEdge4");
460  }
461 
462  //----------------------------------------------------------
463 
464  void RecoJMShower::beginRun(art::Run& /*run*/)
465  {
466  if(!fLoaded){
467  art::ServiceHandle<geo::Geometry> geom;
468  fDetid = geom->DetId();
469  if(pRecoPID==true)LoadDedxHistograms(geom->DetId());
470  LoadTreeWeights(geom->DetId());
471  fLoaded = true;
472  }
473 
474 
475  }
476 
477  void RecoJMShower::postBeginRun(art::Run const& run)
478  {
479  }
480 
481  void RecoJMShower::preEvent(art::Event const& evt) {
482  }
483 
484  //------------------------------------------------------------
485  // Load Probability Histograms and Training Tree Weight Files
486  //------------------------------------------------------------
487 
488  bool RecoJMShower::LoadDedxHistograms(int detid)
489  {
490 
491  std::string libpath = util::EnvExpansion(fLibPath);
492 
493  std::string fnameE(libpath);
494  std::string fnameG(libpath);
495  std::string fnameMu(libpath);
496  std::string fnamePi0(libpath);
497 // std::string fnameHad(libpath);
498  std::string fnameP(libpath);
499  std::string fnameN(libpath);
500  std::string fnamePi(libpath);
501  std::string fnameEqe(libpath);
502  std::string fnameEres(libpath);
503  std::string fnameEdis(libpath);
504  std::string fnameEcoh(libpath);
505  std::string fnameEsg(libpath);
506 
507  if(pPIDHist==1){
508  switch (detid) {
509  case novadaq::cnv::kFARDET:
510 
511  fnameE += "e10kdedxhist_fd_nue_cc.root";
512  fnameG += "g10kdedxhist_fd_nue_cc.root";
513  fnameMu += "mu10kdedxhist_fd_nue_cc.root";
514  fnamePi0 += "pi010kdedxhist_fd_nue_cc.root";
515  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
516  fnameP += "protondedxhist_fd.root";
517  fnameN += "neutrondedxhist_fd.root";
518  fnamePi += "pidedxhist_fd.root";
519  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
520  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
521  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
522  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
523  fnameEsg += "singleehist_fd.root";
524 
525 /*
526  fnameE += "e10kdedxhist_fd_nue_cc.root";
527  fnameG += "g10kdedxhist_genie_fd_nhc.root";
528  fnameMu += "mu10kdedxhist_genie_fd_nhc.root";
529  fnamePi0 += "pi010kdedxhist_genie_fd_nhc.root";
530  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
531  fnameP += "protondedxhist_fd.root";
532  fnameN += "neutrondedxhist_fd.root";
533  fnamePi += "pidedxhist_fd.root";
534  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
535  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
536  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
537  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
538 */
539  break;
540  case novadaq::cnv::kNDOS:
541  fnameE += "e10kdedxhist_fd_nue_cc.root";
542  fnameG += "g10kdedxhist_fd_nue_cc.root";
543  fnameMu += "mu10kdedxhist_fd_nue_cc.root";
544  fnamePi0 += "pi010kdedxhist_fd_nue_cc.root";
545  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
546  fnameP += "protondedxhist_fd.root";
547  fnameN += "neutrondedxhist_fd.root";
548  fnamePi += "pidedxhist_fd.root";
549  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
550  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
551  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
552  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
553  fnameEsg += "singleehist_fd.root";
554  break;
555  case novadaq::cnv::kNEARDET:
556 
557  fnameE += "e10kdedxhist_fd_nue_cc.root";
558  fnameG += "g10kdedxhist_fd_nue_cc.root";
559  fnameMu += "mu10kdedxhist_fd_nue_cc.root";
560  fnamePi0 += "pi010kdedxhist_fd_nue_cc.root";
561  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
562  fnameP += "protondedxhist_fd.root";
563  fnameN += "neutrondedxhist_fd.root";
564  fnamePi += "pidedxhist_fd.root";
565  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
566  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
567  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
568  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
569  fnameEsg += "singleehist_fd.root";
570 
571  break;
572  default:
573  throw cet::exception("BAD DETID") << __FILE__ << ":" << __LINE__
574  << " Bad detector id = " << detid;
575  }
576  }else {
577  switch (detid) {
578  case novadaq::cnv::kFARDET:
579 
580  fnameE += "e10kdedxhist_fd_nue_cc.root";
581  fnameG += "g10kdedxhist_fd_nue_cc.root";
582  fnameMu += "mu10kdedxhist_fd_nue_cc.root";
583  fnamePi0 += "pi010kdedxhist_fd_nue_cc.root";
584  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
585  fnameP += "protondedxhist_fd.root";
586  fnameN += "neutrondedxhist_fd.root";
587  fnamePi += "pidedxhist_fd.root";
588  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
589  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
590  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
591  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
592  fnameEsg += "singleehist_fd.root";
593 
594  break;
595  case novadaq::cnv::kNDOS:
596 
597  fnameE += "e10kdedxhist_fd_nue_cc.root";
598  fnameG += "g10kdedxhist_fd_nue_cc.root";
599  fnameMu += "mu10kdedxhist_fd_nue_cc.root";
600  fnamePi0 += "pi010kdedxhist_fd_nue_cc.root";
601  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
602  fnameP += "protondedxhist_fd.root";
603  fnameN += "neutrondedxhist_fd.root";
604  fnamePi += "pidedxhist_fd.root";
605  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
606  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
607  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
608  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
609  fnameEsg += "singleehist_fd.root";
610 
611  break;
612  case novadaq::cnv::kNEARDET:
613 
614  fnameE += "e10kdedxhist_fd_nue_cc.root";
615  fnameG += "g10kdedxhist_fd_nue_cc.root";
616  fnameMu += "mu10kdedxhist_fd_nue_cc.root";
617  fnamePi0 += "pi010kdedxhist_fd_nue_cc.root";
618  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
619  fnameP += "protondedxhist_fd.root";
620  fnameN += "neutrondedxhist_fd.root";
621  fnamePi += "pidedxhist_fd.root";
622  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
623  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
624  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
625  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
626  fnameEsg += "singleehist_fd.root";
627  break;
628  default:
629  throw cet::exception("BAD DETID") << __FILE__ << ":" << __LINE__
630  << " Bad detector id = " << detid;
631  }
632  }
633 
634 /* histograms stored in /nova/data/users/bianjm/share/development/RecoJMShower/histograms for the test versions
635 
636  std::string fFnameE;
637  std::string fFnameG;
638  std::string fFnameMu;
639  std::string fFnamePi0;
640  std::string fFnameHad;
641  std::string fFnameP;
642  std::string fFnameN;
643  std::string fFnamePi;
644  std::string fFnameEqe;
645  std::string fFnameEres;
646  std::string fFnameEdis;
647  std::string fFnameEcoh;
648 // constructor decides if initialized value is a path or an environment variable
649  cet::search_path sp("FW_SEARCH_PATH");
650  sp.find_file(fnameE, fFnameE);
651  sp.find_file(fnameG, fFnameG);
652  sp.find_file(fnameMu, fFnameMu);
653  sp.find_file(fnamePi0, fFnamePi0);
654  sp.find_file(fnameHad, fFnameHad);
655  sp.find_file(fnameP, fFnameP);
656  sp.find_file(fnameN, fFnameN);
657  sp.find_file(fnamePi, fFnamePi);
658  sp.find_file(fnameEqe, fFnameEqe);
659  sp.find_file(fnameEres, fFnameEres);
660  sp.find_file(fnameEdis, fFnameEdis);
661  sp.find_file(fnameEcoh, fFnameEcoh);
662 
663  struct stat sb;
664 
665  if (fFnameE.empty() || stat(fFnameE.c_str(), &sb)!=0){
666  // Failed to resolve the file name
667  mf::LogWarning("RecoJMShower") << "Histogram -- E " << fFnameE << " not found.";
668  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
669  }
670 
671  if (fFnameEqe.empty() || stat(fFnameEqe.c_str(), &sb)!=0){
672  // Failed to resolve the file name
673  mf::LogWarning("RecoJMShower") << "Histogram -- Eqe " << fFnameEqe << " not found.";
674  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
675  }
676  if (fFnameEres.empty() || stat(fFnameEres.c_str(), &sb)!=0){
677  // Failed to resolve the file name
678  mf::LogWarning("RecoJMShower") << "Histogram -- Eres " << fFnameEres << " not found.";
679  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
680  }
681  if (fFnameEdis.empty() || stat(fFnameEdis.c_str(), &sb)!=0){
682  // Failed to resolve the file name
683  mf::LogWarning("RecoJMShower") << "Histogram -- Edis " << fFnameEdis << " not found.";
684  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
685  }
686  if (fFnameEcoh.empty() || stat(fFnameEcoh.c_str(), &sb)!=0){
687  // Failed to resolve the file name
688  mf::LogWarning("RecoJMShower") << "Histogram -- Ecoh " << fFnameEcoh << " not found.";
689  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
690  }
691 
692  if (fFnameG.empty() || stat(fFnameG.c_str(), &sb)!=0){
693  // Failed to resolve the file name
694  mf::LogWarning("RecoJMShower") << "Histogram -- G " << fFnameG << " not found.";
695  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
696  }
697  if (fFnameMu.empty() || stat(fFnameMu.c_str(), &sb)!=0){
698  // Failed to resolve the file name
699  mf::LogWarning("RecoJMShower") << "Histogram -- Mu " << fFnameMu << " not found.";
700  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
701  }
702  if (fFnamePi0.empty() || stat(fFnamePi0.c_str(), &sb)!=0){
703  // Failed to resolve the file name
704  mf::LogWarning("RecoJMShower") << "Histogram -- Pi0 " << fFnamePi0 << " not found.";
705  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
706  }
707 
708  if (fFnameHad.empty() || stat(fFnameHad.c_str(), &sb)!=0){
709  // Failed to resolve the file name
710  mf::LogWarning("RecoJMShower") << "Histogram -- Had " << fFnameHad << " not found.";
711  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
712  }
713 
714  if (fFnameP.empty() || stat(fFnameP.c_str(), &sb)!=0){
715  // Failed to resolve the file name
716  mf::LogWarning("RecoJMShower") << "Histogram -- P " << fFnameP << " not found.";
717  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
718  }
719 
720  if (fFnameN.empty() || stat(fFnameN.c_str(), &sb)!=0){
721  // Failed to resolve the file name
722  mf::LogWarning("RecoJMShower") << "Histogram -- N " << fFnameN << " not found.";
723  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
724  }
725 
726  if (fFnamePi.empty() || stat(fFnamePi.c_str(), &sb)!=0){
727  // Failed to resolve the file name
728  mf::LogWarning("RecoJMShower") << "Histogram -- Pi " << fFnamePi << " not found.";
729  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
730  }
731 
732 */
733  TFile hfileE(fnameE.c_str());
734  TFile hfileG(fnameG.c_str());
735  TFile hfileMu(fnameMu.c_str());
736  TFile hfilePi0(fnamePi0.c_str());
737 // TFile hfileHad(fnameHad.c_str());
738  TFile hfileP(fnameP.c_str());
739  TFile hfileN(fnameN.c_str());
740  TFile hfilePi(fnamePi.c_str());
741  TFile hfileEqe(fnameEqe.c_str());
742  TFile hfileEres(fnameEres.c_str());
743  TFile hfileEdis(fnameEdis.c_str());
744  TFile hfileEcoh(fnameEcoh.c_str());
745  TFile hfileEsg(fnameEsg.c_str());
746  std::cout<<"PID Dedx File "<<fnameE.c_str()<<" is loaded."<<std::endl;
747  std::cout<<"PID Dedx File "<<fnameG.c_str()<<" is loaded."<<std::endl;
748  std::cout<<"PID Dedx File "<<fnameMu.c_str()<<" is loaded."<<std::endl;
749  std::cout<<"PID Dedx File "<<fnamePi0.c_str()<<" is loaded."<<std::endl;
750 // std::cout<<"PID Dedx File "<<fnameHad.c_str()<<" is loaded."<<std::endl;
751  std::cout<<"PID Dedx File "<<fnameP.c_str()<<" is loaded."<<std::endl;
752  std::cout<<"PID Dedx File "<<fnameN.c_str()<<" is loaded."<<std::endl;
753  std::cout<<"PID Dedx File "<<fnamePi.c_str()<<" is loaded."<<std::endl;
754  std::cout<<"PID Dedx File "<<fnameEqe.c_str()<<" is loaded."<<std::endl;
755  std::cout<<"PID Dedx File "<<fnameEres.c_str()<<" is loaded."<<std::endl;
756  std::cout<<"PID Dedx File "<<fnameEdis.c_str()<<" is loaded."<<std::endl;
757  std::cout<<"PID Dedx File "<<fnameEcoh.c_str()<<" is loaded."<<std::endl;
758  std::cout<<"PID Dedx File "<<fnameEsg.c_str()<<" is loaded."<<std::endl;
759 
760  //The detector is divided into 4 regions in XY, these regions extend in z.
761  //Region 1 is the top left quadrant in XY when looking at the front of the detector
762  //The training histograms are then divided into 11 energy regions from 0 to 5 GeV.
763  //In the training each shower is put into a bin based on its reco energy
764  //The first and last energy bin are 0-0.25 and 4.75-5 GeV. All bins inbetween are 0.5 GeV wide.
765  //For each region and each energy there is a de/dx histogram for each plane from 0 to 200
766  //the histograms being loaded below are for longitudinal dedx.
767  //These training histograms are also broken up by particle type (electron, muon, gamma, proton, pi0, pion, neutron)
768 
769  //NOTE: Jianming is exploring using the shower length as a training variable.
770  //In order to do this the variables below are initialized and then store the
771  //shower length at each energy in each region. This is done by finding the
772  //last plane in which there was energy deposited.
773  //These variables are not currently being used in the production version of the algorithm
774  for(int r=0; r<4; r++){
775  for(int i=0; i<11; i++){
776  fHtExpPlaneE[r][i]=0;
777  fHtExpPlaneG[r][i]=0;
778  fHtExpPlaneMu[r][i]=0;
779  fHtExpPlanePi0[r][i]=0;
780  fHtExpPlaneHad[r][i]=0;
781  fHtExpPlaneP[r][i]=0;
782  fHtExpPlaneN[r][i]=0;
783  fHtExpPlanePi[r][i]=0;
784  fHtExpPlaneEqe[r][i]=0;
785  fHtExpPlaneEres[r][i]=0;
786  fHtExpPlaneEdis[r][i]=0;
787  fHtExpPlaneEcoh[r][i]=0;
788  }
789  }
790 
791  for(int r=0; r<4; r++){
792  for(int i=0; i<41; i++){
793  fHtExpPlaneEsg[r][i]=0;
794  }
795  }
796 
797 
798 
799  for(int r=0; r<4; r++){
800  for(int i=0; i<11; i++){
801  for(int l=0; l<200; l++){
802  char cht[200];
803  sprintf(cht,"ht_%d_%d_%d",r,i,l);
804  fHtPlaneDedxE[r][i][l] = (TH1D*)hfileE.Get(cht);
805  fHtPlaneDedxG[r][i][l] = (TH1D*)hfileG.Get(cht);
806  fHtPlaneDedxMu[r][i][l] = (TH1D*)hfileMu.Get(cht);
807  fHtPlaneDedxPi0[r][i][l] = (TH1D*)hfilePi0.Get(cht);
808  fHtPlaneDedxP[r][i][l] = (TH1D*)hfileP.Get(cht);
809  fHtPlaneDedxN[r][i][l] = (TH1D*)hfileN.Get(cht);
810  fHtPlaneDedxPi[r][i][l] = (TH1D*)hfilePi.Get(cht);
811  fHtPlaneDedxEqe[r][i][l] = (TH1D*)hfileEqe.Get(cht);
812  fHtPlaneDedxEres[r][i][l] = (TH1D*)hfileEres.Get(cht);
813  fHtPlaneDedxEdis[r][i][l] = (TH1D*)hfileEdis.Get(cht);
814  fHtPlaneDedxEcoh[r][i][l] = (TH1D*)hfileEcoh.Get(cht);
815 
816 
817  if(fHtPlaneDedxE[r][i][l]->GetEntries()>0&&fHtExpPlaneE[r][i]==0){
818  if(fHtPlaneDedxE[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxE[r][i][l]->GetEntries()>0.997){
819  fHtExpPlaneE[r][i]=l;
820  }
821  }
822 
823  if(fHtPlaneDedxG[r][i][l]->GetEntries()>0&&fHtExpPlaneG[r][i]==0){
824  if(fHtPlaneDedxG[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxG[r][i][l]->GetEntries()>0.997){
825  fHtExpPlaneG[r][i]=l;
826  }
827  }
828 
829  if(fHtPlaneDedxMu[r][i][l]->GetEntries()>0&&fHtExpPlaneMu[r][i]==0){
830  if(fHtPlaneDedxMu[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxMu[r][i][l]->GetEntries()>0.997){
831  fHtExpPlaneMu[r][i]=l;
832  }
833  }
834 
835  if(fHtPlaneDedxPi0[r][i][l]->GetEntries()>0&&fHtExpPlanePi0[r][i]==0){
836  if(fHtPlaneDedxPi0[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxPi0[r][i][l]->GetEntries()>0.997){
837  fHtExpPlanePi0[r][i]=l;
838  }
839  }
840 
841  if(fHtPlaneDedxP[r][i][l]->GetEntries()>0&&fHtExpPlaneP[r][i]==0){
842  if(fHtPlaneDedxP[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxP[r][i][l]->GetEntries()>0.997){
843  fHtExpPlaneP[r][i]=l;
844  }
845  }
846 
847  if(fHtPlaneDedxN[r][i][l]->GetEntries()>0&&fHtExpPlaneN[r][i]==0){
848  if(fHtPlaneDedxN[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxN[r][i][l]->GetEntries()>0.997){
849  fHtExpPlaneN[r][i]=l;
850  }
851  }
852 
853  if(fHtPlaneDedxPi[r][i][l]->GetEntries()>0&&fHtExpPlanePi[r][i]==0){
854  if(fHtPlaneDedxPi[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxPi[r][i][l]->GetEntries()>0.997){
855  fHtExpPlanePi[r][i]=l;
856  }
857  }
858 
859  if(fHtPlaneDedxEqe[r][i][l]->GetEntries()>0&&fHtExpPlaneEqe[r][i]==0){
860  if(fHtPlaneDedxEqe[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxEqe[r][i][l]->GetEntries()>0.997){
861  fHtExpPlaneEqe[r][i]=l;
862  }
863  }
864 
865  if(fHtPlaneDedxEres[r][i][l]->GetEntries()>0&&fHtExpPlaneEres[r][i]==0){
866  if(fHtPlaneDedxEres[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxEres[r][i][l]->GetEntries()>0.997){
867  fHtExpPlaneEres[r][i]=l;
868  }
869  }
870 
871  if(fHtPlaneDedxEdis[r][i][l]->GetEntries()>0&&fHtExpPlaneEdis[r][i]==0){
872  if(fHtPlaneDedxEdis[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxEdis[r][i][l]->GetEntries()>0.997){
873  fHtExpPlaneEdis[r][i]=l;
874  }
875  }
876 
877  if(fHtPlaneDedxEcoh[r][i][l]->GetEntries()>0&&fHtExpPlaneEcoh[r][i]==0){
878  if(fHtPlaneDedxEcoh[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxEcoh[r][i][l]->GetEntries()>0.997){
879  fHtExpPlaneEcoh[r][i]=l;
880  }
881  }
882  }
883 
884 
885 /*
886  if(r==0){
887  std::cout<<"igev "<<i<<" "<<std::endl;
888  std::cout<<"exp plane E "<<fHtExpPlaneE[r][i]<<std::endl;
889  std::cout<<"exp plane G "<<fHtExpPlaneG[r][i]<<std::endl;
890  std::cout<<"exp plane Mu "<<fHtExpPlaneMu[r][i]<<std::endl;
891  std::cout<<"exp plane Pi0 "<<fHtExpPlanePi0[r][i]<<std::endl;
892  std::cout<<"exp plane Had "<<fHtExpPlaneHad[r][i]<<std::endl;
893  std::cout<<"exp plane P "<<fHtExpPlaneP[r][i]<<std::endl;
894  std::cout<<"exp plane N "<<fHtExpPlaneN[r][i]<<std::endl;
895  std::cout<<"exp plane Pi "<<fHtExpPlanePi[r][i]<<std::endl;
896  std::cout<<"exp plane Eqe "<<fHtExpPlaneEqe[r][i]<<std::endl;
897  std::cout<<"exp plane Eres "<<fHtExpPlaneEres[r][i]<<std::endl;
898  std::cout<<"exp plane Edis "<<fHtExpPlaneEdis[r][i]<<std::endl;
899  std::cout<<"exp plane Ecoh "<<fHtExpPlaneEcoh[r][i]<<std::endl;
900  std::cout<<""<<std::endl;
901  }
902 */
903  }
904  }
905 
906  for(int r=0; r<4; r++){
907  for(int i=0; i<41; i++){
908  for(int l=0; l<200; l++){
909  char cht[200];
910  sprintf(cht,"ht_%d_%d_%d",r,i,l);
911  fHtPlaneDedxEsg[r][i][l] = (TH1D*)hfileEsg.Get(cht);
912  if(fHtPlaneDedxEsg[r][i][l]->GetEntries()>0&&fHtExpPlaneEsg[r][i]==0){
913  if(fHtPlaneDedxEsg[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxEsg[r][i][l]->GetEntries()>0.997){
914  fHtExpPlaneEsg[r][i]=l;
915  }
916  }
917  }
918  }
919  }
920 
921 
922 
923 
924 /*
925  for(int i=0; i<11; i++){
926  for(int l=0; l<200; l++){
927  char cht[200];
928  sprintf(cht,"ht_%d_%d",i,l);
929  fHtPlaneDedxE[i][l] = (TH1D*)hfileE.Get(cht);
930  fHtPlaneDedxG[i][l] = (TH1D*)hfileG.Get(cht);
931  fHtPlaneDedxMu[i][l] = (TH1D*)hfileMu.Get(cht);
932  fHtPlaneDedxPi0[i][l] = (TH1D*)hfilePi0.Get(cht);
933 // fHtPlaneDedxHad[i][l] = (TH1D*)hfileHad.Get(cht);
934  fHtPlaneDedxP[i][l] = (TH1D*)hfileP.Get(cht);
935  fHtPlaneDedxN[i][l] = (TH1D*)hfileN.Get(cht);
936  fHtPlaneDedxPi[i][l] = (TH1D*)hfilePi.Get(cht);
937  fHtPlaneDedxEqe[i][l] = (TH1D*)hfileEqe.Get(cht);
938  fHtPlaneDedxEres[i][l] = (TH1D*)hfileEres.Get(cht);
939  fHtPlaneDedxEdis[i][l] = (TH1D*)hfileEdis.Get(cht);
940  fHtPlaneDedxEcoh[i][l] = (TH1D*)hfileEcoh.Get(cht);
941 
942  char chtmip[200];
943  sprintf(chtmip,"htmip_%d_%d",i,l);
944  fHtmipPlaneDedxE[i][l] = (TH1D*)hfileE.Get(chtmip);
945  fHtmipPlaneDedxG[i][l] = (TH1D*)hfileG.Get(chtmip);
946  fHtmipPlaneDedxMu[i][l] = (TH1D*)hfileMu.Get(chtmip);
947  fHtmipPlaneDedxPi0[i][l] = (TH1D*)hfilePi0.Get(chtmip);
948 // fHtmipPlaneDedxHad[i][l] = (TH1D*)hfileHad.Get(chtmip);
949  fHtmipPlaneDedxP[i][l] = (TH1D*)hfileP.Get(chtmip);
950  fHtmipPlaneDedxN[i][l] = (TH1D*)hfileN.Get(chtmip);
951  fHtmipPlaneDedxPi[i][l] = (TH1D*)hfilePi.Get(chtmip);
952  fHtmipPlaneDedxEqe[i][l] = (TH1D*)hfileEqe.Get(chtmip);
953  fHtmipPlaneDedxEres[i][l] = (TH1D*)hfileEres.Get(chtmip);
954  fHtmipPlaneDedxEdis[i][l] = (TH1D*)hfileEdis.Get(chtmip);
955  fHtmipPlaneDedxEcoh[i][l] = (TH1D*)hfileEcoh.Get(chtmip);
956 
957  char chtshower[200];
958  sprintf(chtshower,"htshower_%d_%d",i,l);
959  fHtshowerPlaneDedxE[i][l] = (TH1D*)hfileE.Get(chtshower);
960  fHtshowerPlaneDedxG[i][l] = (TH1D*)hfileG.Get(chtshower);
961  fHtshowerPlaneDedxMu[i][l] = (TH1D*)hfileMu.Get(chtshower);
962  fHtshowerPlaneDedxPi0[i][l] = (TH1D*)hfilePi0.Get(chtshower);
963 // fHtshowerPlaneDedxHad[i][l] = (TH1D*)hfileHad.Get(chtshower);
964  fHtshowerPlaneDedxP[i][l] = (TH1D*)hfileP.Get(chtshower);
965  fHtshowerPlaneDedxN[i][l] = (TH1D*)hfileN.Get(chtshower);
966  fHtshowerPlaneDedxPi[i][l] = (TH1D*)hfilePi.Get(chtshower);
967  fHtshowerPlaneDedxEqe[i][l] = (TH1D*)hfileEqe.Get(chtshower);
968  fHtshowerPlaneDedxEres[i][l] = (TH1D*)hfileEres.Get(chtshower);
969  fHtshowerPlaneDedxEdis[i][l] = (TH1D*)hfileEdis.Get(chtshower);
970  fHtshowerPlaneDedxEcoh[i][l] = (TH1D*)hfileEcoh.Get(chtshower);
971  }
972  }
973 */
974 
975  //now the histrograms for transverse dedx are also being loaded.
976  for(int r=0; r<4; r++){
977  for(int i=0; i<11; i++){
978  for(int l=0; l<20; l++){
979  char cht[200];
980  sprintf(cht,"httrans_%d_%d_%d",r,i,l);
981  fHtCellTransDedxE[r][i][l] = (TH1D*)hfileE.Get(cht);
982  fHtCellTransDedxG[r][i][l] = (TH1D*)hfileG.Get(cht);
983  fHtCellTransDedxMu[r][i][l] = (TH1D*)hfileMu.Get(cht);
984  fHtCellTransDedxPi0[r][i][l] = (TH1D*)hfilePi0.Get(cht);
985 // fHtCellTransDedxHad[r][i][l] = (TH1D*)hfileHad.Get(cht);
986  fHtCellTransDedxP[r][i][l] = (TH1D*)hfileP.Get(cht);
987  fHtCellTransDedxN[r][i][l] = (TH1D*)hfileN.Get(cht);
988  fHtCellTransDedxPi[r][i][l] = (TH1D*)hfilePi.Get(cht);
989  fHtCellTransDedxEqe[r][i][l] = (TH1D*)hfileEqe.Get(cht);
990  fHtCellTransDedxEres[r][i][l] = (TH1D*)hfileEres.Get(cht);
991  fHtCellTransDedxEdis[r][i][l] = (TH1D*)hfileEdis.Get(cht);
992  fHtCellTransDedxEcoh[r][i][l] = (TH1D*)hfileEcoh.Get(cht);
993  }
994  }
995  }
996 
997  for(int r=0; r<4; r++){
998  for(int i=0; i<41; i++){
999  for(int l=0; l<20; l++){
1000  char cht[200];
1001  sprintf(cht,"httrans_%d_%d_%d",r,i,l);
1002  fHtCellTransDedxEsg[r][i][l] = (TH1D*)hfileEsg.Get(cht);
1003  }
1004  }
1005  }
1006 
1007 
1008  hfileE.Close();
1009  hfileG.Close();
1010  hfileMu.Close();
1011  hfilePi0.Close();
1012 // hfileHad.Close();
1013  hfileP.Close();
1014  hfileN.Close();
1015  hfilePi.Close();
1016  hfileEqe.Close();
1017  hfileEres.Close();
1018  hfileEdis.Close();
1019  hfileEcoh.Close();
1020  hfileEsg.Close();
1021  return true;
1022  }
1023 
1024  bool RecoJMShower::LoadTreeWeights(int detid)
1025  {
1026  std::string libpath = util::EnvExpansion(fLibPath);
1027  std::string fnameWeightsShape(libpath);
1028  std::string fnameWeightsShapeElec(libpath);
1029  std::string fnameWeightsShapeE(libpath);
1030  std::string fnameWeightsShapeNoCos(libpath);
1031  std::string fnameWeightsShapeENoCos(libpath);
1032  std::string fnameWeightsShapeEL(libpath);
1033  std::string fnameWeightsShapeepi0(libpath);
1034  std::string fnameWeightsShapeepi0EL(libpath);
1035 
1036 
1037 // std::string fnameWeightsShapeQE(libpath);
1038 // std::string fnameWeightsShapeRES(libpath);
1039 // std::string fnameWeightsShapeDIS(libpath);
1040  switch (detid) {
1041  case novadaq::cnv::kFARDET:
1042  fnameWeightsShape += "weights_shape_fd_elec.txt";
1043  fnameWeightsShapeElec += "weights_shower_fd_elec.txt";
1044  fnameWeightsShapeE += "weights_shape_fd_elec_E.txt";
1045  fnameWeightsShapeNoCos += "weights_shape_fd_elec_NoCos.txt";
1046  fnameWeightsShapeENoCos += "weights_shape_fd_elec_E_NoCos.txt";
1047  fnameWeightsShapeEL += "weights_shape_fd_elec_E_EL.txt";
1048  fnameWeightsShapeepi0 += "weights_shape_fd_elec_E_epi0.txt";
1049  fnameWeightsShapeepi0EL += "weights_shape_fd_elec_E_epi0EL.txt";
1050 // fnameWeightsShapeQE += "weights_shape_fd_shower_QE.txt";
1051 // fnameWeightsShapeRES += "weights_shape_fd_shower_RES.txt";
1052 // fnameWeightsShapeDIS += "weights_shape_fd_shower_DIS.txt";
1053  break;
1054  case novadaq::cnv::kNDOS:
1055  fnameWeightsShape += "weights_shape_fd_elec.txt";
1056  fnameWeightsShapeElec += "weights_shower_fd_elec.txt";
1057  fnameWeightsShapeE += "weights_shape_fd_elec_E.txt";
1058  fnameWeightsShapeNoCos += "weights_shape_fd_elec_NoCos.txt";
1059  fnameWeightsShapeENoCos += "weights_shape_fd_elec_E_NoCos.txt";
1060  fnameWeightsShapeEL += "weights_shape_fd_elec_E_EL.txt";
1061  fnameWeightsShapeepi0 += "weights_shape_fd_elec_E_epi0.txt";
1062  fnameWeightsShapeepi0EL += "weights_shape_fd_elec_E_epi0EL.txt";
1063 // fnameWeightsShapeQE += "weights_shape_fd_shower_QE.txt";
1064 // fnameWeightsShapeRES += "weights_shape_fd_shower_RES.txt";
1065 // fnameWeightsShapeDIS += "weights_shape_fd_shower_DIS.txt";
1066  break;
1067  case novadaq::cnv::kNEARDET:
1068  fnameWeightsShape += "weights_shape_fd_elec.txt";
1069  fnameWeightsShapeElec += "weights_shower_fd_elec.txt";
1070  fnameWeightsShapeE += "weights_shape_fd_elec_E.txt";
1071  fnameWeightsShapeNoCos += "weights_shape_fd_elec_NoCos.txt";
1072  fnameWeightsShapeENoCos += "weights_shape_fd_elec_E_NoCos.txt";
1073  fnameWeightsShapeEL += "weights_shape_fd_elec_E_EL.txt";
1074  fnameWeightsShapeepi0 += "weights_shape_fd_elec_E_epi0.txt";
1075  fnameWeightsShapeepi0EL += "weights_shape_fd_elec_E_epi0EL.txt";
1076 // fnameWeightsShapeQE += "weights_shape_fd_shower_QE.txt";
1077 // fnameWeightsShapeRES += "weights_shape_fd_shower_RES.txt";
1078 // fnameWeightsShapeDIS += "weights_shape_fd_shower_DIS.txt";
1079  break;
1080  default:
1081  throw cet::exception("BAD DETID") << __FILE__ << ":" << __LINE__
1082  << " Bad detector id = " << detid;
1083  }
1084 /* Histograms stored in /nova/data/users/bianjm/share/development/RecoJMShower/histograms for this test versions
1085 
1086  cet::search_path sp("FW_SEARCH_PATH");
1087  sp.find_file(fnameWeightsShape, fFnameWeightsShape);
1088  struct stat sb;
1089  if (fFnameWeightsShape.empty() || stat(fFnameWeightsShape.c_str(), &sb)!=0){
1090  // Failed to resolve the file name
1091  mf::LogWarning("RecoJMShower") << "Training Weight -- " << fFnameWeightsShape << " not found.";
1092  throw cet::exception("FILE_NOT_FOUND") << "ANN Weight file not found.";
1093  }
1094 */
1095 // std::cout<<"Trained ANN "<<fnameWeightsShape.c_str()<<" is loaded."<<std::endl;
1096 
1097  double egLLL;
1098  double egLLT;
1099  double emuLLL;
1100  double emuLLT;
1101  double epi0LLL;
1102  double epi0LLT;
1103  double epLLL;
1104  double epLLT;
1105  double enLLL;
1106  double enLLT;
1107  double epiLLL;
1108  double epiLLT;
1109  double dedx0;
1110  double dedx1;
1111  double dedx2;
1112  double dedx3;
1113  double dedx4;
1114  double dedx5;
1115 
1116  double vtxgev;
1117  double vtxcentergev;
1118  double vtxcentergev5cell;
1119  double vtxcentergev10cell;
1120  double vtxcentergev15cell;
1121  double vtxcentergev20cell;
1122  double pi0mass;
1123  double shE;
1124  double nuE;
1125  double gap;
1126  double length;
1127  double costheta;
1128  int type;
1129 
1130  TTree *trainTree_shape = new TTree("trainTree_shape","m_trainTree_shape");
1131 
1132  trainTree_shape->Branch("egLLL", &egLLL, "egLLL/D");
1133  trainTree_shape->Branch("egLLT", &egLLT, "egLLT/D");
1134  trainTree_shape->Branch("emuLLL", &emuLLL, "emuLLL/D");
1135  trainTree_shape->Branch("emuLLT", &emuLLT, "emuLLT/D");
1136  trainTree_shape->Branch("epi0LLL", &epi0LLL, "epi0LLL/D");
1137  trainTree_shape->Branch("epi0LLT", &epi0LLT, "epi0LLT/D");
1138  trainTree_shape->Branch("epLLL", &epLLL, "epLLL/D");
1139  trainTree_shape->Branch("epLLT", &epLLT, "epLLT/D");
1140  trainTree_shape->Branch("enLLL", &enLLL, "enLLL/D");
1141  trainTree_shape->Branch("enLLT", &enLLT, "enLLT/D");
1142  trainTree_shape->Branch("epiLLL", &epiLLL, "epiLLL/D");
1143  trainTree_shape->Branch("epiLLT", &epiLLT, "epiLLT/D");
1144 
1145  trainTree_shape->Branch("dedx0", &dedx0, "dedx0/D");
1146  trainTree_shape->Branch("dedx1", &dedx1, "dedx1/D");
1147  trainTree_shape->Branch("dedx2", &dedx2, "dedx2/D");
1148  trainTree_shape->Branch("dedx3", &dedx3, "dedx3/D");
1149  trainTree_shape->Branch("dedx4", &dedx4, "dedx4/D");
1150  trainTree_shape->Branch("dedx5", &dedx5, "dedx5/D");
1151 
1152 /*
1153  trainTree_shape->Branch("eelgLLL", &eelgLLL, "eelgLLL/D");
1154  trainTree_shape->Branch("eelgLLT", &eelgLLT, "eelgLLT/D");
1155  trainTree_shape->Branch("eelmuLLL", &eelmuLLL, "eelmuLLL/D");
1156  trainTree_shape->Branch("eelmuLLT", &eelmuLLT, "eelmuLLT/D");
1157  trainTree_shape->Branch("eelpi0LLL", &eelpi0LLL, "eelpi0LLL/D");
1158  trainTree_shape->Branch("eelpi0LLT", &eelpi0LLT, "eelpi0LLT/D");
1159  trainTree_shape->Branch("eelpLLL", &eelpLLL, "eelpLLL/D");
1160  trainTree_shape->Branch("eelpLLT", &eelpLLT, "eelpLLT/D");
1161  trainTree_shape->Branch("eelnLLL", &eelnLLL, "eelnLLL/D");
1162  trainTree_shape->Branch("eelnLLT", &eelnLLT, "eelnLLT/D");
1163  trainTree_shape->Branch("eelpiLLL", &eelpiLLL, "eelpiLLL/D");
1164  trainTree_shape->Branch("eelpiLLT", &eelpiLLT, "eelpiLLT/D");
1165 */
1166 
1167  trainTree_shape->Branch("gap", &gap, "gap/D");
1168  trainTree_shape->Branch("pi0mass", &pi0mass, "pi0mass/D");
1169  trainTree_shape->Branch("vtxgev", &vtxgev, "vtxgev/D");
1170  trainTree_shape->Branch("shE", &shE, "shE/D");
1171  trainTree_shape->Branch("nuE", &nuE, "nuE/D");
1172  trainTree_shape->Branch("costheta", &costheta, "costheta/D");
1173  trainTree_shape->Branch("length", &length, "length/D");
1174  trainTree_shape->Branch("type", &type, "type/I");
1175 
1176  trainTree_shape->Branch("vtxcentergev", &vtxcentergev, "vtxcentergev/D");
1177  trainTree_shape->Branch("vtxcentergev5cell", &vtxcentergev5cell, "vtxcentergev5cell/D");
1178  trainTree_shape->Branch("vtxcentergev10cell", &vtxcentergev10cell, "vtxcentergev10cell/D");
1179  trainTree_shape->Branch("vtxcentergev15cell", &vtxcentergev15cell, "vtxcentergev15cell/D");
1180  trainTree_shape->Branch("vtxcentergev20cell", &vtxcentergev20cell, "vtxcentergev20cell/D");
1181 
1182 
1183  mlp_shape = new TMultiLayerPerceptron("@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@pi0mass,@shE,@vtxgev,@gap,@costheta:22:12:6:type", trainTree_shape);
1184  mlp_shape->LoadWeights(fnameWeightsShape.c_str());
1185 
1186  mlp_shape_NoCos = new TMultiLayerPerceptron("@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@pi0mass,@shE,@vtxgev,@gap:22:12:6:type", trainTree_shape);
1187  mlp_shape_NoCos->LoadWeights(fnameWeightsShapeNoCos.c_str());
1188 
1189 
1190  mlp_shape_elec = new TMultiLayerPerceptron("@egLLL,@egLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@shE,@gap:18:12:6:type", trainTree_shape);
1191 
1192  mlp_shape_elec->LoadWeights(fnameWeightsShapeElec.c_str());
1193 
1194  mlp_shape_E = new TMultiLayerPerceptron("@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@pi0mass,@shE,@vtxgev,@gap,@costheta,@nuE:22:12:6:type", trainTree_shape);
1195  mlp_shape_E->LoadWeights(fnameWeightsShapeE.c_str());
1196 
1197 
1198  mlp_shape_E_NoCos = new TMultiLayerPerceptron("@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@pi0mass,@shE,@vtxgev,@gap,@nuE:22:12:6:type", trainTree_shape);
1199  mlp_shape_E_NoCos->LoadWeights(fnameWeightsShapeENoCos.c_str());
1200 
1201 
1202  mlp_shape_EL = new TMultiLayerPerceptron("@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT:22:12:6:type", trainTree_shape);
1203  mlp_shape_EL->LoadWeights(fnameWeightsShapeEL.c_str());
1204 
1205 
1206  mlp_shape_epi0 = new TMultiLayerPerceptron("@dedx0,@dedx1,@dedx2,@dedx3,@costheta:25:12:6:type", trainTree_shape);
1207  mlp_shape_epi0->LoadWeights(fnameWeightsShapeepi0.c_str());
1208 
1209 
1210  mlp_shape_epi0EL = new TMultiLayerPerceptron("@dedx0,@dedx1,@dedx2,@dedx3:25:12:6:type", trainTree_shape);
1211  mlp_shape_epi0EL->LoadWeights(fnameWeightsShapeepi0EL.c_str());
1212 
1213 /*
1214  mlp_shape = new TMultiLayerPerceptron("@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@pi0mass,@shE,@gap,@costheta,@vtxgev,@length:25:16:9:type", trainTree_shape);
1215  mlp_shape->LoadWeights(fnameWeightsShape.c_str());
1216 
1217 
1218  mlp_shape = new TMultiLayerPerceptron("egLLL,egLLT,epi0LLL,epi0LLT,epLLL,epLLT,enLLL,enLLT,epiLLL,epiLLT,pi0mass,vtxgev:18:12:6:type",trainTree_shape);
1219  mlp_shape->LoadWeights(fnameWeightsShape.c_str());
1220 */
1221 
1222 /*
1223  mlp_shape_QE = new TMultiLayerPerceptron("@egLLL,@egLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@gap:18:12:6:type",trainTree_shape);
1224  mlp_shape_RES = new TMultiLayerPerceptron("@egLLL,@egLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@gap:18:12:6:type",trainTree_shape);
1225  mlp_shape_DIS = new TMultiLayerPerceptron("@egLLL,@egLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@gap:18:12:6:type",trainTree_shape);
1226  mlp_shape_QE->LoadWeights(fnameWeightsShapeQE.c_str());
1227  mlp_shape_RES->LoadWeights(fnameWeightsShapeRES.c_str());
1228  mlp_shape_DIS->LoadWeights(fnameWeightsShapeDIS.c_str());
1229 */
1230  return true;
1231  }
1232 
1233  void RecoJMShower::produce(art::Event& evt)
1234  {
1235  isrealdata = evt.isRealData();
1236  fprongidx.clear();
1237 
1238  //containers to store shower objects in event and the associations to slices
1239  std::unique_ptr< std::vector<jmshower::JMShower> > recoshowercol(new std::vector<jmshower::JMShower> );
1240  std::unique_ptr< std::vector<jmshower::JMShower> > evtrecoshowercol(new std::vector<jmshower::JMShower> );
1241  std::unique_ptr< art::Assns<jmshower::JMShower, rb::Cluster> > assns(new art::Assns<jmshower::JMShower, rb::Cluster>);
1242  std::unique_ptr< art::Assns<jmshower::JMShower, rb::Prong> > prongassns(new art::Assns<jmshower::JMShower, rb::Prong>);
1243  // TODO: Prong-jmshower association is temporary. It should be deleted once ShowerLID validation is done.
1244 
1245  //pass the event to the function that does the work, creates a vector of shower objects
1246  //pass the event to the function that does the work, creates a vector of shower objects
1247 // std::vector<jmshower::JMShower> jmshowercol = RecoJMShower::RecoShowers(evt);
1248  std::vector<jmshower::JMShower> jmshowercol = RecoJMShower::RecoShowers(*this, evt, recoshowercol, prongassns);
1249 
1250  //get the slices contained in each event
1251  art::Handle<std::vector<rb::Cluster> > slicecol;
1252  evt.getByLabel(pSliceDir, slicecol);
1253  //loop over the event showers and put them in the event
1254  for(unsigned int i=0; i<jmshowercol.size();i++){
1255  evtrecoshowercol->push_back(jmshowercol[i]);
1256  //each shower has stored the slice it comes from, pull this out now to make proper association
1257  unsigned int iSlice = jmshowercol[i].ISlice();
1258  util::CreateAssn(*this,evt,*evtrecoshowercol,art::Ptr<rb::Cluster>(slicecol, iSlice),*assns);
1259  }
1260 
1261  evt.put(std::move(evtrecoshowercol));
1262  evt.put(std::move(assns));
1263  evt.put(std::move(prongassns));
1264  return;
1265  }
1266 
1267 
1268 // std::vector<jmshower::JMShower> RecoJMShower::RecoShowers(const art::Event& evt){
1269  std::vector<jmshower::JMShower> RecoJMShower::RecoShowers(art::EDProducer const& prod,
1270  art::Event& evt,
1271  std::unique_ptr< std::vector<jmshower::JMShower> >&recoshowercol,
1272  std::unique_ptr< art::Assns<jmshower::JMShower, rb::Prong> >& prongassns){
1273 
1274 // int run = evt.run();
1275 // int srun = evt.subRun();
1276 // int event = evt.id().event();
1277  //load services
1278  art::ServiceHandle<geo::Geometry> geom;
1279  art::ServiceHandle<calib::Calibrator> cal;
1280 
1281 //Event level
1282  art::Handle<std::vector<rb::Cluster> > slicecol;
1283  evt.getByLabel(pSliceDir, slicecol);
1284 
1285  std::vector<jmshower::JMShower> evtshowercol;
1286  evtshowercol.clear();
1287 
1288  //get verticies associated with each slice
1289  art::FindManyP<rb::Vertex> fmv(slicecol, evt, pVertexLabel); // get vertices
1290 //std::cout<<"pTrackDir , pInstLabel "<<pTrackDir<<" "<<pInstLabel<<std::endl;
1291 // art::FindManyP<rb::Prong> fmprong(slicecol, evt, pTrackDir); // get prong
1292  art::FindManyP<rb::Prong> fmprong(slicecol, evt, art::InputTag(pTrackDir, pInstLabel));
1293  if(pRecoVtx==true)if (!(fmv.isValid())) return evtshowercol;
1294  if (!(fmprong.isValid())) return evtshowercol;
1295  for(unsigned int ic = 0; ic < slicecol->size(); ++ic){
1296  if(rb::IsFiltered(evt,slicecol,ic)) continue;
1297  art::Ptr<rb::Cluster> clust(slicecol,ic);
1298  if(clust->IsNoise()) continue;
1299  std::vector<jmshower::JMShower> showercol;
1300  showercol.clear();
1301 
1302  std::vector<art::Ptr<rb::Prong> > prong2d3d = fmprong.at(ic);
1303  std::vector<art::Ptr<rb::Prong> > track;
1304  track.clear();
1305  for(unsigned int iPro = 0; iPro < prong2d3d.size(); ++iPro){
1306  if( prong2d3d[iPro]->Is3D())track.push_back(prong2d3d[iPro]);
1307  fprongidx.push_back(iPro);
1308  }
1309  //if there are no 3D prongs or crazy # of prongs skip this slice.
1310  //NOTE: This is the only form of Preselection in RecoJMShower, it must
1311  //have a 3D track to be able to proceed.
1312 //std::cout<<"n fuz "<<track.size()<<std::endl;
1313  if (track.size()==0||track.size()>9999){
1314  continue;
1315  }
1316 
1317  //map of cell hits to the slice they came from
1318  //map of cell hits to track they came from
1319  std::map<geo::OfflineChan, std::vector<int>> hittrkmap;
1320  std::map<geo::OfflineChan, std::vector<int>> hitshmap;
1321 
1322  //containers to slice and track info
1323  art::PtrVector< rb::CellHit > cellhitlist;
1324  std::vector<const rb::Cluster* > cluster;
1325  std::vector<rb::Prong> vCluster;
1326  std::vector<rb::Prong> vXCluster;
1327  std::vector<rb::Prong> vYCluster;
1328  std::vector<bool> vNantag;// A tage to deal with crazy prongs
1329  std::vector<double> vtxGeV;
1330  std::vector<double> vtxCenterGeV;
1331  std::vector<double> trkGeV;
1332  std::vector<double> vtxCenterGeV5cell;
1333  std::vector<double> vtxCenterGeV10cell;
1334  std::vector<double> vtxCenterGeV15cell;
1335  std::vector<double> vtxCenterGeV20cell;
1336 
1337  double sliceGeV = 0;
1338  double sliceEtot = 0;
1339  double sliceEdge10GeV = 0;//slice energy close to detector edge - 10 cell/plane
1340  double sliceEdge5GeV = 0;//slice energy close to detector edge - 5 cell/plane
1341  double sliceEdge2GeV = 0;//slice energy close to detector edge - 2 cell/plane
1342 
1343  std::vector<int> shNCellToEdge;//Mininum distance from shower cells to detector edges
1344  std::vector<int> shClosetEdge;//Closet edge for a shower
1345  std::vector<double> shEdge20GeV;//shower energy close to detector edge - 20 cell/plane
1346  std::vector<double> shEdge15GeV;//shower energy close to detector edge - 15 cell/plane
1347  std::vector<double> shEdge10GeV;//shower energy close to detector edge - 10 cell/plane
1348  std::vector<double> shEdge5GeV;//shower energy close to detector edge - 5 cell/plane
1349  std::vector<double> shEdge2GeV;//shower energy close to detector edge - 2 cell/plane
1350 
1351  double sliceExclGeV = 0;//slice energy close to detector edge - 2 cell/plane
1352 
1353  rb::Cluster excludeCluster; //cluster of cells not in any tracks in slice
1354 
1355  //clear containers
1356  hittrkmap.clear();
1357  hitshmap.clear();
1358  fHitDistMap.clear();
1359  fHitEWeightMap.clear();
1360  fHitEMap.clear();
1361 
1362  fShowerPClusterCol.clear(); //longitudinal shower cluster
1363  fTrackPClusterCol.clear(); //longitudinal track cluster
1364  fShowerTCClusterCol.clear(); //transverse shower cluster
1365  fExcludeShowerCol.clear(); //shower objects for hits in slice not put into a track
1366 
1367  //clear containers
1368  cellhitlist.clear();
1369  cluster.clear();
1370  vCluster.clear();
1371  vXCluster.clear();
1372  vYCluster.clear();
1373  vNantag.clear();
1374  vtxGeV.clear();
1375  trkGeV.clear();
1376  shNCellToEdge.clear();
1377  shClosetEdge.clear();
1378  shEdge20GeV.clear();
1379  shEdge15GeV.clear();
1380  shEdge10GeV.clear();
1381  shEdge5GeV.clear();
1382  shEdge2GeV.clear();
1383 
1384  excludeCluster.Clear();
1385 
1386  //get the slices
1387 
1388  std::vector<TVector3> evtvtxcol; //hold vevent vertex
1389  std::vector<unsigned int> evtvtxplanecol; //plane associated with vertex
1390 
1391  //clear containers
1392  evtvtxcol.clear();
1393  evtvtxplanecol.clear();
1394 
1395  //initialize vectors to hold vertex info for each view
1396  Double_t cellXYZ1[3]={0,0,0};
1397  Double_t cellXYZ2[3]={0,0,0};
1398  geom->Plane(0)->Cell(0)->GetCenter(cellXYZ1);
1399  geom->Plane(1)->Cell(0)->GetCenter(cellXYZ2);
1400 
1401  //if this option is set, find the elastic arms vertex for each slice
1402  if(pRecoVtx==true){
1403  //NOTE: if there is no vertex, then there are no tracks, cannot do EID, so return
1404  //This is a form of "preselection"
1405  //if there was nue preselection, skip slice if it failed
1406  std::vector<art::Ptr<rb::Vertex> > v = fmv.at(ic);
1407  //loop over verticies associated with slice
1408  for (unsigned int j=0; j<v.size(); ++j) {
1409  TVector3 evtVtx(0,0,0);
1410  unsigned int vtxPlane = 0;
1411  //store vertex coordinates
1412  evtVtx.SetX(v[j]->GetX());
1413  evtVtx.SetY(v[j]->GetY());
1414  evtVtx.SetZ(v[j]->GetZ());
1415 
1416  //find plane number closest to vertex
1417  geo::CellUniqueId cid;
1418 
1419  for(int trial = 0; trial < 4; ++trial){
1420  // Figure out what cell the vertex was in. Step 3cm in z in case of
1421  // failure. If that didn't work, it must be in plastic between cells,
1422  // move laterally in x or y.
1423 
1424  double offsetX = 0, offsetY = 0;
1425  if(trial == 1) offsetX = -2;
1426  if(trial == 2) offsetY = -2;
1427  if(trial == 3) offsetX = offsetY = -2;
1428 
1429  cid = geom->CellId(v[j]->GetX()+offsetX, v[j]->GetY()+offsetY, v[j]->GetZ(),
1430  0, 0, 1, 3);
1431  if(cid)
1432  break;
1433  }
1434 
1435  if( cid )
1436  vtxPlane = geom->getPlaneID( cid );
1437  else
1438  vtxPlane = clust->MinPlane();
1439 
1440  //store vertex, plane of vertex, and slice corresponding to vertex
1441  evtvtxcol.push_back(evtVtx);
1442  evtvtxplanecol.push_back(vtxPlane);
1443  }
1444  }
1445 
1446  //Clear hittrkmap
1447  for(unsigned int i=0;i<clust->NXCell();i++){
1448  geo::OfflineChan key(clust->XCell(i)->Plane(),clust->XCell(i)->Cell());
1449  hittrkmap[key].clear(); //set match between hit and track to default
1450  hitshmap[key].clear(); //set match between hit and shower to default
1451  }
1452 
1453  //NOTE: When using FuzzyKVertex tracks cell hits can belong to more then one track
1454  //as this is written now the hit will be associated with the last track - has been solved by using a map
1455  for(unsigned int i=0;i<track.size();i++){ //loop over tracks
1456  for(unsigned int ii=0;ii<track[i]->NCell();ii++){ //loop over cells in tracks
1457  geo::OfflineChan key(track[i]->Cell(ii)->Plane(),track[i]->Cell(ii)->Cell()); //make key value
1458  hittrkmap[key].push_back((int)i); //store track associated with cell
1459  hitshmap[key].push_back((int)i); //store shower associated with cell
1460  }
1461  }
1462 
1463  //TODO: This is setting the transverse coordinate of a view to the coordinate of the
1464  //first cell in that view. This would be better if the average transverse coordinate was used.
1465  //This needs to be corrected but would also require retraining, needs to be looked at.
1466  const double wx = (clust->NXCell() > 0) ? clust->W(clust->XCell(0).get()) : 0;
1467  const double wy = (clust->NYCell() > 0) ? clust->W(clust->YCell(0).get()) : 0;
1468 
1469  //loop over x cells
1470  for(unsigned int i=0;i<clust->NXCell();i++){
1471  geo::OfflineChan key(clust->XCell(i)->Plane(),clust->XCell(i)->Cell());
1472  //only use cell hits on a track or cell hits with pe above threshold level set if fcl file
1473  if((clust->XCell(i)->PE())>pCellHitPHThresh||hittrkmap[key].size()>0){
1474  //make reco hits out of cells to sum energy
1475  const art::Ptr<rb::CellHit>& chit = clust->XCell(i);
1476  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, wx ));
1477 
1478  double gev=0;
1479  geo::View_t view = geom->Plane(chit->Plane())->View();
1480 
1481  double pigtail = geom->GetPigtail(geom->Plane(chit->Plane())->Cell(chit->Cell())->Id());
1482  //double wxtoelec = wx - (pigtail - 92.6624);//W - (pigtail - avg. pigtail)
1483  //double wytoelec = wy - (pigtail - 92.6624);//W - (pigtail - avg. pigtail)
1484 
1485  double w1 = wx;
1486  if (view == geo::kY) w1 = wy;
1487  double xyzd1[4];
1488  calib::GetXYZD(chit->OfflineChan(), w1, xyzd1);
1489  double readDist = xyzd1[3];
1490 
1491  if(rhit.IsCalibrated()){
1492  gev = rhit.GeV();
1493  }else if(pUseUncalib == true){// use MC calibration if a cell is not calibrated
1494 // std::cout<<"slice, plane, cell, wx, gev "<<ic<<" "<<chit->Plane()<<" "<<chit->Cell()<<" "<<wx<<" Not calibrated?!" << std::endl;
1495  double adc = double(chit->ADC());
1496  double Att = 1;
1497  if(view == geo::kX){
1498  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
1499  }else if(view == geo::kY){
1500  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
1501  }
1502  gev = adc/Att;
1503  }
1504 
1505  double ecor = 1.;
1506  double disttoelec = geom->Plane(chit->Plane())->Cell(chit->Cell())->HalfL()+pigtail-wx;
1507  if(evt.isRealData()==true&&pCorDataCell)ecor = RecoJMShower::CellEDataMC(gev, disttoelec, view, fDetid); // correct Data/MC difference for cell energy
1508  gev = gev/ecor;
1509 
1510  sliceGeV += gev;
1511  if(std::abs((int)chit->Cell()-(int)pXEdgeId1)<10 ||
1512  std::abs((int)chit->Cell()-(int)pXEdgeId2)<10 ||
1513  std::abs((int)chit->Cell()-(int)pXEdgeId3)<10 ||
1514  std::abs((int)chit->Cell()-(int)pXEdgeId4)<10 ||
1515  std::abs((int)chit->Plane()-(int)pZEdgeId1)<10 ||
1516  std::abs((int)chit->Plane()-(int)pZEdgeId2)<10 ||
1517  std::abs((int)chit->Plane()-(int)pZEdgeId3)<10 ||
1518  std::abs((int)chit->Plane()-(int)pZEdgeId4)<10){
1519  sliceEdge10GeV += gev;
1520  }
1521 
1522  if(std::abs((int)chit->Cell()-(int)pXEdgeId1)<5 ||
1523  std::abs((int)chit->Cell()-(int)pXEdgeId2)<5 ||
1524  std::abs((int)chit->Cell()-(int)pXEdgeId3)<5 ||
1525  std::abs((int)chit->Cell()-(int)pXEdgeId4)<5 ||
1526  std::abs((int)chit->Plane()-(int)pZEdgeId1)<5 ||
1527  std::abs((int)chit->Plane()-(int)pZEdgeId2)<5 ||
1528  std::abs((int)chit->Plane()-(int)pZEdgeId3)<5 ||
1529  std::abs((int)chit->Plane()-(int)pZEdgeId4)<5){
1530  sliceEdge5GeV += gev;
1531  }
1532 
1533  if(std::abs((int)chit->Cell()-(int)pXEdgeId1)<2 ||
1534  std::abs((int)chit->Cell()-(int)pXEdgeId2)<2 ||
1535  std::abs((int)chit->Cell()-(int)pXEdgeId3)<2 ||
1536  std::abs((int)chit->Cell()-(int)pXEdgeId4)<2 ||
1537  std::abs((int)chit->Plane()-(int)pZEdgeId1)<2 ||
1538  std::abs((int)chit->Plane()-(int)pZEdgeId2)<2 ||
1539  std::abs((int)chit->Plane()-(int)pZEdgeId3)<2 ||
1540  std::abs((int)chit->Plane()-(int)pZEdgeId4)<2){
1541  sliceEdge2GeV += gev;
1542  }
1543  cellhitlist.push_back(clust->XCell(i));
1544  }
1545  }
1546 
1547  //loop over y cells
1548  for(unsigned int i=0;i<clust->NYCell();i++){
1549  geo::OfflineChan key(clust->YCell(i)->Plane(),clust->YCell(i)->Cell());
1550  //use hits above pe threshold or on a track
1551  if((clust->YCell(i)->PE())>pCellHitPHThresh||hittrkmap[key].size()>0){
1552  //make reco hits out of cells
1553  const art::Ptr<rb::CellHit>& chit = clust->YCell(i);
1554  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, wy ));
1555  double gev=0;
1556  geo::View_t view = geom->Plane(chit->Plane())->View();
1557  double pigtail = geom->GetPigtail(geom->Plane(chit->Plane())->Cell(chit->Cell())->Id());
1558  //double wxtoelec = wx - (pigtail - 92.6624);//W - (pigtail - avg. pigtail)
1559  //double wytoelec = wy - (pigtail - 92.6624);//W - (pigtail - avg. pigtail)
1560 
1561  double w1 = wx;
1562  if (view == geo::kY) w1 = wy;
1563  double xyzd1[4];
1564  calib::GetXYZD(chit->OfflineChan(), w1, xyzd1);
1565  double readDist = xyzd1[3];
1566 
1567  if(rhit.IsCalibrated()){
1568  gev = rhit.GeV();
1569  }else if(pUseUncalib == true){// use MC calibration if a cell is not calibrated
1570  double adc = double(chit->ADC());
1571  double Att = 1;
1572  if(view == geo::kX){
1573  Att = RecoJMShower::CellADCToGeV(readDist, view, novadaq::cnv::kFARDET);
1574  }else if(view == geo::kY){
1575  Att = RecoJMShower::CellADCToGeV(readDist, view, novadaq::cnv::kFARDET);
1576  }
1577  gev = adc/Att;
1578  }
1579 
1580  double ecor = 1.;
1581  double disttoelec = geom->Plane(chit->Plane())->Cell(chit->Cell())->HalfL()+pigtail-wy;
1582  if(evt.isRealData()==true&&pCorDataCell)ecor = RecoJMShower::CellEDataMC(gev, disttoelec, view, fDetid); // correct Data/MC difference for cell energy
1583  gev = gev/ecor;
1584 
1585  sliceGeV += gev;
1586 
1587  if(std::abs((int)chit->Cell()-(int)pYEdgeId1)<10 ||
1588  std::abs((int)chit->Cell()-(int)pYEdgeId2)<10 ||
1589  std::abs((int)chit->Cell()-(int)pYEdgeId3)<10 ||
1590  std::abs((int)chit->Cell()-(int)pYEdgeId4)<10 ||
1591  std::abs((int)chit->Plane()-(int)pZEdgeId1)<10 ||
1592  std::abs((int)chit->Plane()-(int)pZEdgeId2)<10 ||
1593  std::abs((int)chit->Plane()-(int)pZEdgeId3)<10 ||
1594  std::abs((int)chit->Plane()-(int)pZEdgeId4)<10){
1595  sliceEdge10GeV += gev;
1596  }
1597  if(std::abs((int)chit->Cell()-(int)pYEdgeId1)<5 ||
1598  std::abs((int)chit->Cell()-(int)pYEdgeId2)<5 ||
1599  std::abs((int)chit->Cell()-(int)pYEdgeId3)<5 ||
1600  std::abs((int)chit->Cell()-(int)pYEdgeId4)<5 ||
1601  std::abs((int)chit->Plane()-(int)pZEdgeId1)<5 ||
1602  std::abs((int)chit->Plane()-(int)pZEdgeId2)<5 ||
1603  std::abs((int)chit->Plane()-(int)pZEdgeId3)<5 ||
1604  std::abs((int)chit->Plane()-(int)pZEdgeId4)<5){
1605  sliceEdge5GeV += gev;
1606  }
1607  if(std::abs((int)chit->Cell()-(int)pYEdgeId1)<2 ||
1608  std::abs((int)chit->Cell()-(int)pYEdgeId2)<2 ||
1609  std::abs((int)chit->Cell()-(int)pYEdgeId3)<2 ||
1610  std::abs((int)chit->Cell()-(int)pYEdgeId4)<2 ||
1611  std::abs((int)chit->Plane()-(int)pZEdgeId1)<2 ||
1612  std::abs((int)chit->Plane()-(int)pZEdgeId2)<2 ||
1613  std::abs((int)chit->Plane()-(int)pZEdgeId3)<2 ||
1614  std::abs((int)chit->Plane()-(int)pZEdgeId4)<2){
1615  sliceEdge2GeV += gev;
1616  }
1617 
1618  cellhitlist.push_back(clust->YCell(i));
1619  }
1620  }
1621 
1622 
1623  //declare and clear some containers for shower variables
1624  std::vector<TVector3> startcol,endcol,dircol;
1625  std::vector<TVector3> shstartcol,shendcol;
1626  std::vector<unsigned int> startplanecol,stopplanecol;
1627  std::vector<unsigned int> shstartplanecol,shstopplanecol;
1628  std::vector<bool> backtagcol;
1629  std::vector<double> gapcol;
1630  std::vector<double> vtxdocacol;
1631  startcol.clear();
1632  endcol.clear();
1633  dircol.clear();
1634  startplanecol.clear();
1635  stopplanecol.clear();
1636  backtagcol.clear();
1637  gapcol.clear();
1638  vtxdocacol.clear();
1639  shstartcol.clear();
1640  shendcol.clear();
1641  shstartplanecol.clear();
1642  shstopplanecol.clear();
1643 
1644 
1645  //loop over tracks
1646  //this loop stores infomation about the start and stop of each track
1647  //and the distance to the vertex
1648  for(unsigned int i=0;i<track.size();i++){
1649  //initialize cell container
1650  art::PtrVector<rb::CellHit> trackcells;
1651  trackcells.clear();
1652  //copy the track into a new track object
1653  rb::Prong Clust(trackcells,track[i]->Start(),track[i]->Dir());
1654  // Clust.AppendTrajectoryPoint(GetShwStop(track[i]));
1655 
1656  //push this copy into vector that can be manipulated
1657 
1658  TVector3 tStart = track[i]->Start();
1659  bool nantag = false;// reset track start/stop if nan vale is found
1660  if(std::isnan(tStart.x())||std::abs(tStart.x())>10000){tStart.SetX(track[i]->MinX());nantag = true;}
1661  if(std::isnan(tStart.y())||std::abs(tStart.y())>10000){tStart.SetY(track[i]->MinY());nantag = true;}
1662  if(std::isnan(tStart.z())||std::abs(tStart.z())>10000){tStart.SetZ(track[i]->MinZ());nantag = true;}
1663  TVector3 tStop = RecoJMShower::GetShwStop(track[i]);
1664  if(std::isnan(tStop.x()||std::abs(tStop.x())>10000)){tStop.SetX(track[i]->MaxX());nantag = true;}
1665  if(std::isnan(tStop.y()||std::abs(tStop.y())>10000)){tStop.SetY(track[i]->MaxY());nantag = true;}
1666  if(std::isnan(tStop.z()||std::abs(tStop.z())>10000)){tStop.SetZ(track[i]->MaxZ());nantag = true;}
1667  TVector3 tVect = (tStop-tStart).Unit();
1668  if(nantag == true){
1669  // Clust.ClearTrajectoryPoints();
1670  Clust.SetStart(tStart);
1671  // Clust.AppendTrajectoryPoint(tStop);
1672  Clust.SetDir(tVect);
1673  }
1674 
1675  vCluster.push_back(Clust);
1676  vXCluster.push_back(Clust);
1677  vYCluster.push_back(Clust);
1678  vNantag.push_back(nantag);
1679 
1680  vtxGeV.push_back(0); //initialize vertex energy
1681  vtxCenterGeV.push_back(0); //initialize vertex energy
1682  trkGeV.push_back(0); //initialize track energy
1683  shNCellToEdge.push_back(9999);//initialize mininum distance from shower cells to detector edges
1684  shClosetEdge.push_back(-1);//initialize Closet edge for a shower
1685  shEdge20GeV.push_back(0);//initialize shower energy close to detector edge - 20 cell/plane
1686  shEdge15GeV.push_back(0);//initialize shower energy close to detector edge - 15 cell/plane
1687  shEdge10GeV.push_back(0);//initialize shower energy close to detector edge - 10 cell/plane
1688  shEdge5GeV.push_back(0);//initialize shower energy close to detector edge - 5 cell/plane
1689  shEdge2GeV.push_back(0);//initialize shower energy close to detector edge - 2 cell/plane
1690 
1691  vtxCenterGeV5cell.push_back(0); //initialize vertex energy
1692  vtxCenterGeV10cell.push_back(0); //initialize vertex energy
1693  vtxCenterGeV15cell.push_back(0); //initialize vertex energy
1694  vtxCenterGeV20cell.push_back(0); //initialize vertex energy
1695 
1696  //initialize event vertex to start of track and first plane in track
1697  TVector3 evtVtx(track[i]->Start()[0],track[i]->Start()[1],track[i]->Start()[2]);
1698  geo::OfflineChan key(track[i]->Cell(0)->Plane(),track[i]->Cell(0)->Cell());
1699  //if the vertex is being used by the EID (Which is the default option and what the training is based on)
1700  //then find the vertex in the slice the track comes from
1701  if(pRecoVtx==true&&pVtxFit==true){
1702  for(unsigned iv=0;iv<evtvtxcol.size();iv++){
1703  // double doca = RecoJMShower::GetPointDoca(evtvtxcol[iv], track[i]->Start(), RecoJMShower::GetShwStop(track[i]));
1704  // std::cout<<"dist,doca "<<(evtVtx-evtvtxcol[iv]).Mag()<<" "<<doca<<" "<<std::endl;
1705  // if((evtVtx-evtvtxcol[iv]).Mag()<60.&&doca<20.){
1706  evtVtx = evtvtxcol[iv]; //get the vertex for that slice
1707  break;
1708  // }
1709  }
1710  }
1711 
1712  //get the distance of closest approach to vertex
1713  double vtxDoca = RecoJMShower::GetPointDoca(evtVtx, track[i]->Start(), RecoJMShower::GetShwStop(track[i]));
1714  double vtxToStart = (evtVtx-track[i]->Start()).Mag();
1715  double vtxToStop = (evtVtx-RecoJMShower::GetShwStop(track[i])).Mag();
1716  double vtxIntToStart = sqrt(vtxToStart*vtxToStart-vtxDoca*vtxDoca);
1717  double vtxIntToStop = sqrt(vtxToStop*vtxToStop-vtxDoca*vtxDoca);
1718  TVector3 vtxStart = track[i]->Start();
1719  TVector3 vtxStop = RecoJMShower::GetShwStop(track[i]);
1720  // unsigned int vtxStopPlane = track[i]->MaxPlane();
1721 
1722  //record track start and end planes
1723  TVector3 shStart = track[i]->Start();
1724  TVector3 shStop = RecoJMShower::GetShwStop(track[i]);
1725 
1726  unsigned int shStartPlane = track[i]->MinPlane();
1727  unsigned int shStopPlane = track[i]->MaxPlane();
1728  double vtxTrkLength = (vtxStop-vtxStart).Mag();
1729  bool showerVtxBackTag = false;
1730  //flip start and stop for backward tracks if needed
1731  if(vtxToStart>vtxToStop) {
1732  showerVtxBackTag = true;
1733  vtxStart = RecoJMShower::GetShwStop(track[i]);
1734  vtxStop = track[i]->Start();
1735  double tmp = vtxIntToStart;
1736  vtxIntToStart = vtxIntToStop;
1737  vtxIntToStop = tmp;
1738  shStart = RecoJMShower::GetShwStop(track[i]);
1739  shStop = track[i]->Start();
1740  }
1741  //if(std::abs((int)track[i]->MinPlane()-(int)vtxPlane)>std::abs((int)track[i]->MaxPlane()-(int)vtxPlane)){
1742  if(track[i]->Dir()[2] < 0.0){
1743  // vtxStopPlane = track[i]->MinPlane();
1744  shStartPlane = track[i]->MaxPlane();
1745  shStopPlane = track[i]->MinPlane();
1746  }
1747  if(vtxIntToStop<vtxTrkLength){
1748  vtxIntToStart = -vtxIntToStart;
1749  }
1750 
1751  TVector3 trkVtxVect = (vtxStop-evtVtx).Unit();
1752  TVector3 shVect = (shStop-shStart).Unit();
1753 
1754  //Set start and direction
1755  if(pVtxFit==true){
1756  bool vtxnantag = false;// reset track start/stop if nan vale is found
1757  if(std::isnan(tStart.x())||std::abs(tStart.x())>10000){tStart.SetX(track[i]->MinX());vtxnantag = true;}
1758  if(std::isnan(tStart.y())||std::abs(tStart.y())>10000){tStart.SetY(track[i]->MinY());vtxnantag = true;}
1759  if(std::isnan(tStart.z())||std::abs(tStart.z())>10000){tStart.SetZ(track[i]->MinZ());vtxnantag = true;}
1760  if(std::isnan(tStop.x()||std::abs(tStop.x())>10000)){tStop.SetX(track[i]->MaxX());vtxnantag = true;}
1761  if(std::isnan(tStop.y()||std::abs(tStop.y())>10000)){tStop.SetY(track[i]->MaxY());vtxnantag = true;}
1762  if(std::isnan(tStop.z()||std::abs(tStop.z())>10000)){tStop.SetZ(track[i]->MaxZ());vtxnantag = true;}
1763  if(vtxnantag==true) shVect = (shStop-shStart).Unit();
1764  vCluster[i].SetStart(shStart);
1765  vCluster[i].SetDir(shVect);
1766  vXCluster[i].SetStart(shStart);
1767  vXCluster[i].SetDir(shVect);
1768  vYCluster[i].SetStart(shStart);
1769  vYCluster[i].SetDir(shVect);
1770 
1771  }
1772  // reverse shower using vertex
1773  //put track start and end points into vectors
1774  startcol.push_back(shStart);
1775  endcol.push_back(shStop);
1776  dircol.push_back(shVect);
1777  startplanecol.push_back(shStartPlane);
1778  stopplanecol.push_back(shStopPlane);
1779 
1780  backtagcol.push_back(showerVtxBackTag);
1781  //record the gap from vertex to the start of the track
1782  //if the vertex is inside a track record distance of closest approach
1783  //TODO: For fuzzyK there can be an internal gap, where the start of the track is misleading
1784  //Ruth has looked into this, will try to incorporate internal gaps into this algorithm
1785  double shGap = vtxIntToStart;
1786  if(vtxIntToStart>=0&&vtxIntToStart<vtxDoca)shGap = vtxDoca;
1787  if(vtxIntToStart<0)shGap = vtxDoca;
1788 
1789  //push gap and track start/stop variables into vectors
1790  gapcol.push_back(shGap);
1791 // vtxdocacol.push_back(vtxDoca);
1792  shstartcol.push_back(shStart);
1793  shendcol.push_back(shStop);
1794  shstartplanecol.push_back(shStartPlane);
1795  shstopplanecol.push_back(shStopPlane);
1796  }
1797 
1798  //*********************************************
1799  // Cell energy deconvolution and calibration
1800  //*********************************************
1801 
1802  //now looping over all the cell hits in slice to do energy calibration
1803  for (unsigned int i=0;i<cellhitlist.size();i++) {
1804  //get iterator for ith cell
1805  std::vector<art::Ptr<rb::CellHit> >::iterator itr = cellhitlist.begin() + i;
1806  //turn cell into an offlinechan so it can be found in the maps
1807  geo::OfflineChan key((*itr)->Plane(),(*itr)->Cell());
1808  // double cellL = 2*geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->HalfL();// longest
1809  // double cellD = 2*geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->HalfD();// longitudinal
1810  //store the cell width
1811  double cellW = 2*geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->HalfW();
1812  // std::cout<<"cell L,D,W "<<cellL<<" "<<cellD<<" "<<cellW<<std::endl;
1813  //initialize weights
1814  double thiteweighttotal = 0;
1815  double thiteweighttotalatt = 0;
1816  double thiteweight[100];
1817  // double trkattcal[20];
1818 
1819  //now loop over tracks to calculate vertex energy
1820  for(unsigned int itrk=0;itrk<track.size();itrk++){
1821  if(vNantag[itrk]==true)continue;
1822  geo::View_t view = geom->Plane((*itr)->Plane())->View();
1823  bool trktag=false;// on orignal tracks ?
1824  for(unsigned int imap=0;imap<hittrkmap[key].size();imap++){
1825  if(hittrkmap[key][imap]==(int)itrk)trktag=true;
1826  }
1827  if(trktag==true){//store x, y view clusters for track quality check.
1828  if(view == geo::kX)vXCluster[itrk].Add((*itr));
1829  if(view == geo::kY)vYCluster[itrk].Add((*itr));
1830  }
1831 
1832 
1833  //TVector3 trkdir=track[itrk]->Dir();
1834  //TVector3 trkstart=track[itrk]->Start();
1835  //TVector3 trkstop=RecoJMShower::GetShwStop(track[itrk]);
1836 
1837  //get the track start, stop and direction
1838  TVector3 trkdir=dircol[itrk];
1839  TVector3 trkstart=startcol[itrk];
1840  TVector3 trkstop=endcol[itrk];
1841 
1842  //clear that tracks hit weight
1843  thiteweight[itrk]=0;
1844  //trkattcal[itrk] = 0;
1845 
1846  //get the difference between the plane of the cell hit and the first plane of the track
1847  int plane = (int)(*itr)->Plane()-(int)startplanecol[itrk];
1848  if(startplanecol[itrk]>stopplanecol[itrk])plane = -plane; //reverse if track is backwards
1849  if(itrk==0){
1850  for(unsigned int icc=0;icc<track[itrk]->NCell();icc++){
1851  geo::OfflineChan tt(track[itrk]->Cell(icc)->Plane(),track[itrk]->Cell(icc)->Cell());
1852  }
1853  }
1854  //vertex region is 8 planes in front and back of the start of a track
1855  //first looking if cell hit is in backwards region
1856  double celltrkdist= RecoJMShower::GetCellDistToTrk((*itr)->Plane(),(*itr)->Cell(),shstartcol[itrk],shendcol[itrk]);
1857 
1858  if(plane>=-8&&plane<0){
1859  double cellgev = 0;
1860  geo::View_t view = geom->Plane((*itr)->Plane())->View();
1861  TVector3 trkhitpos = RecoJMShower::GetTrkHitPos((*itr)->Plane(),(*itr)->Cell(),track[itrk]);
1862  double w1 = track[itrk]->W((*itr).get());
1863  double xyzd1[4];
1864  calib::GetXYZD((*itr)->OfflineChan(), w1, xyzd1);
1865  double readDist = xyzd1[3];
1866  //compute an attenuation correction based on hit distance to readout. This is done for the FD
1867  //Correction factor determined in MC by comparing simulated calibration to truth
1868 
1869  if(fDetid!=novadaq::cnv::kFARDET||pCorAbsCell==false){
1870  // for NDOS or ND
1871  double witrk = track[itrk]->W((*itr).get());
1872  //To consider extra fiber of each cell in MC, for data dont apply this because cell are calibrated one by one.
1873  double pigtail = geom->GetPigtail(geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->Id());
1874  if(evt.isRealData()==false&&pCorPigtail==true)witrk = witrk - (pigtail - 92.6624);//W - (pigtail - avg. pigtail)
1875  const rb::RecoHit rhit = cal->MakeRecoHit(*(*itr), witrk);
1876  if(rhit.IsCalibrated()){
1877  cellgev = rhit.GeV();
1878  }else if(pUseUncalib == true){
1879  double adc = double((*itr)->ADC());
1880  double Att = 1;
1881  if(view == geo::kX){
1882  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
1883  }else if(view == geo::kY){
1884  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
1885  }
1886  cellgev = adc/Att;
1887  }
1888  }
1889  else{
1890  double adc = double((*itr)->ADC());
1891  double Att = 1;
1892  //Factor is different for each detector view
1893  if(view == geo::kX){
1894  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
1895  }else if(view == geo::kY){
1896  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
1897  }
1898  cellgev = adc/Att;
1899  }
1900  double ecor = 1.0;
1901  if(evt.isRealData()==true&&pCorDataCell)ecor = RecoJMShower::CellEDataMC(cellgev, readDist, view, fDetid); // correct Data/MC difference for cell energy
1902  vtxGeV[itrk]+=cellgev/ecor;
1903  if(celltrkdist<19.5*cellW)vtxCenterGeV[itrk]+=cellgev/ecor;
1904  if(celltrkdist<=5.0*cellW)vtxCenterGeV5cell[itrk]+=cellgev/ecor;
1905  if(celltrkdist<=10.0*cellW)vtxCenterGeV10cell[itrk]+=cellgev/ecor;
1906  if(celltrkdist<=15.0*cellW)vtxCenterGeV15cell[itrk]+=cellgev/ecor;
1907  if(celltrkdist<=20.0*cellW)vtxCenterGeV20cell[itrk]+=cellgev/ecor;
1908  }
1909 
1910  int shnplane = std::abs((int)stopplanecol[i]-(int)startplanecol[i])+1;
1911 
1912  if(plane>=-1*shnplane&&plane<2*shnplane&&celltrkdist<19.5*cellW){// Calculate minimum distance from shower cells to detector edges, extrapolating the shower to -L from start point and +L from stop point. Cells within 20 cells to the core are considered.
1913  double cellgev = 0;
1914  geo::View_t view = geom->Plane((*itr)->Plane())->View();
1915  TVector3 trkhitpos = RecoJMShower::GetTrkHitPos((*itr)->Plane(),(*itr)->Cell(),track[itrk]);
1916  //compute an attenuation correction based on hit distance to readout. This is done for the FD
1917  //Correction factor determined in MC by comparing simulated calibration to truth
1918  double w1 = track[itrk]->W((*itr).get());
1919  double xyzd1[4];
1920  calib::GetXYZD((*itr)->OfflineChan(), w1, xyzd1);
1921  double readDist = xyzd1[3];
1922  if(fDetid!=novadaq::cnv::kFARDET||pCorAbsCell==false){
1923  double witrk = track[itrk]->W((*itr).get());
1924  //To consider extra fiber of each cell in MC, for data dont apply this because cell are calibrated one by one.
1925  double pigtail = geom->GetPigtail(geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->Id());
1926  if(evt.isRealData()==false&&pCorPigtail==true)witrk = witrk - (pigtail - 92.6624);//W - (pigtail - avg. pigtail)
1927  const rb::RecoHit rhit = cal->MakeRecoHit(*(*itr), witrk);
1928  if(rhit.IsCalibrated()){
1929  cellgev = rhit.GeV();
1930  }else if(pUseUncalib == true){
1931  double adc = double((*itr)->ADC());
1932  double Att = 1;
1933  if(view == geo::kX){
1934  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
1935  }else if(view == geo::kY){
1936  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
1937  }
1938  cellgev = adc/Att;
1939  }
1940  }
1941  else{
1942  double adc = double((*itr)->ADC());
1943  double Att = 1;
1944  //Factor is different for each detector view
1945  if(view == geo::kX){
1946  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
1947  }else if(view == geo::kY){
1948  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
1949  }
1950  cellgev = adc/Att;
1951  }
1952  double ecor = 1.0;
1953  if(evt.isRealData()==true&&pCorDataCell)ecor = RecoJMShower::CellEDataMC(cellgev, readDist, view, fDetid); // correct Data/MC difference for cell energy
1954 
1955  int dx1 = std::abs((int)(*itr)->Cell()-(int)pXEdgeId1);
1956  int dx2 = std::abs((int)(*itr)->Cell()-(int)pXEdgeId2);
1957  int dy1 = std::abs((int)(*itr)->Cell()-(int)pYEdgeId1);
1958  int dy2 = std::abs((int)(*itr)->Cell()-(int)pYEdgeId2);
1959  int dz1 = std::abs((int)(*itr)->Plane()-(int)pZEdgeId1);
1960  int dz2 = std::abs((int)(*itr)->Plane()-(int)pZEdgeId2);
1961  int cmd = 9999;
1962  int ce = -1;
1963  if(cmd>dx1){cmd=dx1;ce=1;}
1964  if(cmd>dx2){cmd=dx2;ce=2;}
1965  if(cmd>dy1){cmd=dy1;ce=3;}
1966  if(cmd>dy2){cmd=dy2;ce=4;}
1967  if(cmd>dz1){cmd=dz1;ce=5;}
1968  if(cmd>dz2){cmd=dz2;ce=6;}
1969  if(shNCellToEdge[itrk]>cmd){
1970  shNCellToEdge[itrk]=cmd;
1971  shClosetEdge[itrk]=ce;
1972  }
1973  if(cmd<2)shEdge2GeV[itrk] += cellgev/ecor;
1974  if(cmd<5)shEdge5GeV[itrk] += cellgev/ecor;
1975  if(cmd<10)shEdge10GeV[itrk] += cellgev/ecor;
1976  if(cmd<15)shEdge15GeV[itrk] += cellgev/ecor;
1977  if(cmd<20)shEdge20GeV[itrk] += cellgev/ecor;
1978  }
1979 
1980 
1981  //if the hit occurs before(after) the start(end) of the track continue
1982  if((*itr)->Plane()<std::min(startplanecol[itrk],stopplanecol[itrk])||(*itr)->Plane()>std::max(startplanecol[itrk],stopplanecol[itrk])){fHitDistMap[key].push_back(0);continue;}//for vtx fit
1983  //get the position of hit with respect to track. (distance to readout,transverse coord, transverse dist from hit to track)
1984  TVector3 trkhitpos = RecoJMShower::GetTrkHitPos((*itr)->Plane(),(*itr)->Cell(),track[itrk]);
1985 
1986  double w1 = track[itrk]->W((*itr).get());
1987  double xyzd1[4];
1988  calib::GetXYZD((*itr)->OfflineChan(), w1, xyzd1);
1989  double readDist = xyzd1[3];
1990  // double celltrkdist= trkhitpos.z();
1991  // double celltrkdist= RecoJMShower::GetCellDistToTrk((*itr)->Plane(),(*itr)->Cell(),track[itrk].Start(),RecoJMShower::GetShwStop(track[itrk]));
1992  //get the perpendicular distance from the cell to the track
1993 // double celltrkdist= RecoJMShower::GetCellDistToTrk((*itr)->Plane(),(*itr)->Cell(),shstartcol[itrk],shendcol[itrk]);//orig
1994  /*
1995  if((startplanecol[itrk]<std::min(shstartplanecol[itrk],shstopplanecol[itrk])||startplanecol[itrk]>std::max(shstartplanecol[itrk],shstopplanecol[itrk]))&&((*itr)->Plane()>=std::min(startplanecol[itrk],shstartplanecol[itrk])&&(*itr)->Plane()<=std::max(startplanecol[itrk],shstartplanecol[itrk]))&&(startcol[itrk]-shstartcol[itrk]).Mag()>cellW&&pVtxFit==true){
1996  celltrkdist= RecoJMShower::GetCellDistToTrk((*itr)->Plane(),(*itr)->Cell(),startcol[itrk],shstartcol[itrk]);
1997  }
1998  */
1999 
2000  // double celltrkdistorig = RecoJMShower::GetCellDistToTrk((*itr)->Plane(),(*itr)->Cell(),track[itrk]->Start(),RecoJMShower::GetShwStop(track[itrk]));
2001  //record the perp distance from the cell track inside the map
2002  fHitDistMap[key].push_back(celltrkdist);
2003 
2004  if(plane>=8&&celltrkdist>19.5*cellW&&trktag==false)continue;
2005  //if(plane<8&&celltrkdist>1.5*cellW&&trktag==false){
2006  if(0<=plane&&plane<8&&celltrkdist>1.5*cellW&&trktag==false){
2007  double cellgev = 0;
2008 // geo::View_t view = geom->Plane((*itr)->Plane())->View();
2009 // TVector3 trkhitpos = RecoJMShower::GetTrkHitPos((*itr)->Plane(),(*itr)->Cell(),track[itrk]);
2010  if(fDetid!=novadaq::cnv::kFARDET||pCorAbsCell==false){
2011  // for NDOS or ND
2012  double witrk = track[itrk]->W((*itr).get());
2013  //To consider extra fiber of each cell in MC, for data dont apply this because cell are calibrated one by one.
2014  double pigtail = geom->GetPigtail(geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->Id());
2015  if(evt.isRealData()==false&&pCorPigtail==true)witrk = witrk - (pigtail - 92.6624);//W - (pigtail - avg. pigtail)
2016  const rb::RecoHit rhit = cal->MakeRecoHit(*(*itr), witrk);
2017  if(rhit.IsCalibrated()){
2018  cellgev = rhit.GeV();
2019  }else if(pUseUncalib == true){
2020  double adc = double((*itr)->ADC());
2021  double Att = 1;
2022  if(view == geo::kX){
2023  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2024  }else if(view == geo::kY){
2025  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2026  }
2027  cellgev = adc/Att;
2028  }
2029  }else{
2030  double adc = double((*itr)->ADC());
2031  double Att = 1;
2032  if(view == geo::kX){
2033  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2034  }else if(view == geo::kY){
2035  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2036  }
2037  cellgev = adc/Att;
2038  }
2039  double ecor = 1.0;
2040  if(evt.isRealData()==true&&pCorDataCell)ecor = RecoJMShower::CellEDataMC(cellgev, readDist, view, fDetid); // correct Data/MC difference for cell energy
2041  vtxGeV[itrk]+=cellgev/ecor;
2042  if(celltrkdist<19.5*cellW)vtxCenterGeV[itrk]+=cellgev/ecor;
2043  if(celltrkdist<=5.0*cellW)vtxCenterGeV5cell[itrk]+=cellgev/ecor;
2044  if(celltrkdist<=10.0*cellW)vtxCenterGeV10cell[itrk]+=cellgev/ecor;
2045  if(celltrkdist<=15.0*cellW)vtxCenterGeV15cell[itrk]+=cellgev/ecor;
2046  if(celltrkdist<=20.0*cellW)vtxCenterGeV20cell[itrk]+=cellgev/ecor;
2047 
2048  continue;
2049  }
2050  if(pRecluster==false&&trktag==false)continue;
2051  // if(plane<8&&celltrkdist>19.5*cellW)continue;
2052  vCluster[itrk].Add((*itr));
2053  bool isinsh = false;//Add this cell to cell-shower map
2054  for(unsigned int ish=0;ish<hitshmap[key].size();ish++){
2055  if(hitshmap[key][ish]==(int)itrk)isinsh=true;
2056  }
2057  if(isinsh==false) hitshmap[key].push_back((int)itrk);
2058 
2059  double witrk = track[itrk]->W((*itr).get());
2060  //To conder extra fiber of each cell in MC, for data dont apply this because cell are calibrated one by one.
2061  double pigtail = geom->GetPigtail(geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->Id());
2062  if(evt.isRealData()==false&&pCorPigtail==true)witrk = witrk - (pigtail - 92.6624);//W - (pigtail - avg. pigtail)
2063  double peattitrk = 0;
2064  const rb::RecoHit rhit = cal->MakeRecoHit(*(*itr), witrk);
2065  if(rhit.IsCalibrated()){
2066  peattitrk = rhit.GeV();
2067  }else if(pUseUncalib == true){
2068  double adc = double((*itr)->ADC());
2069  double Att = 1;
2070  if(view == geo::kX){
2071  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2072  }else if(view == geo::kY){
2073  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2074  }
2075  peattitrk = adc/Att;
2076  }
2077 // const double peattitrk = cal->GetPECorr(*(*itr), witrk);
2078  // trkattcal[itrk] = peattitrk;
2079  const double ai = peattitrk/(*itr)->PE();
2080  thiteweight[itrk] = track[itrk]->TotalPE()*exp(-3.28769e-01*celltrkdist);//Shower lateral profile constant from a fit to (Transverse dE/dx)/(shower Energy) as a function of distance to shower core for Monte Carlo electron.
2081  // thiteweight[itrk] = track[itrk]->TotalPE()*exp(-2.1*trkhitpos.z()/8.5);
2082  thiteweighttotal += thiteweight[itrk];
2083  thiteweighttotalatt += thiteweight[itrk]/ai;
2084  }
2085 
2086  for(unsigned int itrk=0;itrk<track.size();itrk++){
2087  if(vNantag[itrk]==true){
2088  fHitEWeightMap[key].push_back(0);
2089  fHitEMap[key]+=0;
2090  continue;
2091  }
2092  // const double gev = trueEnergy[key]*(thiteweight[itrk]/thiteweighttotalatt);
2093  if(thiteweighttotal>0){
2094  double gev = 0;
2095  geo::View_t view = geom->Plane((*itr)->Plane())->View();
2096  TVector3 trkhitpos = RecoJMShower::GetTrkHitPos((*itr)->Plane(),(*itr)->Cell(),track[itrk]);
2097  double w1 = track[itrk]->W((*itr).get());
2098  double xyzd1[4];
2099  calib::GetXYZD((*itr)->OfflineChan(), w1, xyzd1);
2100  double readDist = xyzd1[3];
2101 
2102  // if(fDetid==novadaq::cnv::kFARDET||fDetid==novadaq::cnv::kNDOS||fDetid==novadaq::cnv::kNEARDET)
2103 // double witrk = track[itrk]->W((*itr).get());
2104  if(fDetid!=novadaq::cnv::kFARDET||pCorAbsCell==false){
2105  gev = (*itr)->PE()*(thiteweight[itrk]/thiteweighttotalatt);
2106 // const rb::RecoHit rhit = cal->MakeRecoHit((*itr).get(), witrk);
2107 // if(rhit.IsCalibrated()){
2108 // gev = rhit.GeV();
2109 // }else if(pUseUncalib == true){
2110 // double adc = double((*itr)->ADC());
2111 // double Att = 1;
2112 // if(view == geo::kX){
2113 // Att = RecoJMShower::CellADCToGeV(trkhitpos[0], view, fDetid);
2114 // }else if(view == geo::kY){
2115 // Att = RecoJMShower::CellADCToGeV(trkhitpos[0], view, fDetid);
2116 // }
2117 // gev = adc/Att;
2118 // } // for NDOS or ND
2119  }else{
2120  // for FD, Attenuation effect calibrated for development (Product S12.02.14),
2121  // should not be changed to keep the energy scale identical to the one used in ANN training!
2122  double adc = double((*itr)->ADC());
2123  double Att = 1;
2124  if(view == geo::kX){
2125  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2126  }else if(view == geo::kY){
2127  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2128  }
2129  gev = adc*(thiteweight[itrk]/thiteweighttotal)/Att;
2130  }
2131 
2132  double ecor = 1.0;
2133  if(evt.isRealData()==true&&pCorDataCell)ecor = RecoJMShower::CellEDataMC(gev, readDist, view, fDetid); // correct Data/MC difference for cell energy
2134  gev = gev/ecor;
2135  fHitEWeightMap[key].push_back(gev);
2136  fHitEMap[key]+=gev;
2137  }else{
2138  fHitEWeightMap[key].push_back(0);
2139  fHitEMap[key]+=0;
2140  }
2141  }
2142  }
2143 
2144  for(unsigned int itrk=0;itrk<track.size();itrk++){// calculate orignal prong energy
2145  if(vNantag[itrk]==true)continue;
2146  for(unsigned int i=0;i<track[itrk]->NCell();i++){
2147  art::Ptr< rb::CellHit > trkcell = track[itrk]->Cell(i);
2148  double gev = 0;
2149  geo::View_t view = geom->Plane(trkcell->Plane())->View();
2150  TVector3 trkhitpos = RecoJMShower::GetTrkHitPos(trkcell->Plane(),trkcell->Cell(),track[itrk]);
2151  double w1 = track[itrk]->W(trkcell.get());
2152  double xyzd1[4];
2153  calib::GetXYZD(trkcell->OfflineChan(), w1, xyzd1);
2154  double readDist = xyzd1[3];
2155  if(fDetid==novadaq::cnv::kFARDET||fDetid==novadaq::cnv::kNDOS||fDetid==novadaq::cnv::kNEARDET||pCorAbsCell==false){
2156  // if(fDetid!=novadaq::cnv::kFARDET)
2157  // for NDOS or ND
2158  double witrk = track[itrk]->W(trkcell.get());
2159  //To con sider extra fiber of each cell in MC, for data dont apply this because cell are calibrated one by one.
2160  double pigtail = geom->GetPigtail(geom->Plane(trkcell->Plane())->Cell(trkcell->Cell())->Id());
2161  if(evt.isRealData()==false&&pCorPigtail==true)witrk = witrk - (pigtail - 92.6624);//W - (pigtail - avg. pigtail)
2162  const rb::RecoHit rhit = cal->MakeRecoHit(*trkcell, witrk);
2163  if(rhit.IsCalibrated()){
2164  gev = rhit.GeV();
2165  }else if(pUseUncalib == true){
2166  double adc = double(trkcell->ADC());
2167  double Att = 1;
2168  if(view == geo::kX){
2169  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2170  }else if(view == geo::kY){
2171  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2172  }
2173  gev = adc/Att;
2174  }
2175  }else{
2176  // for FD, Attenuation effect calibrated for development (Product S12.02.14),
2177  // should not be changed to keep the energy scale identical to the one used in ANN training!
2178  double adc = double(trkcell->ADC());
2179  double Att = 1;
2180  if(view == geo::kX){
2181  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2182  }else if(view == geo::kY){
2183  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2184  }
2185  gev = adc/Att;
2186  }
2187  double ecor = 1.0;
2188  if(evt.isRealData()==true&&pCorDataCell)ecor = RecoJMShower::CellEDataMC(gev, readDist, view, fDetid); // correct Data/MC difference for cell energy
2189  gev = gev/ecor;
2190  trkGeV[itrk]+= gev;
2191  }
2192  }
2193 
2194  //loop over x cells
2195  for(unsigned int i=0;i<clust->NXCell();i++){
2196  geo::OfflineChan key(clust->XCell(i)->Plane(),clust->XCell(i)->Cell());
2197  //only use cell hits that not in any showers and with pe above threshold level set if fcl file
2198  if((clust->XCell(i)->PE())>pCellHitPHThresh&&hitshmap[key].size()==0){
2199  //make reco hits out of cells to sum energy
2200  const art::Ptr<rb::CellHit>& chit = clust->XCell(i);
2201  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, wx ));
2202 
2203  double gev=0;
2204  geo::View_t view = geom->Plane(chit->Plane())->View();
2205  double pigtail = geom->GetPigtail(geom->Plane(chit->Plane())->Cell(chit->Cell())->Id());
2206  //double wxtoelec = wx - (pigtail - 92.6624);//W - (pigtail - avg. pigtail)
2207  //double wytoelec = wy - (pigtail - 92.6624);//W - (pigtail - avg. pigtail)
2208  double w1 = wx;
2209  if (view == geo::kY) w1 = wy;
2210  double xyzd1[4];
2211  calib::GetXYZD(chit->OfflineChan(), w1, xyzd1);
2212  double readDist = xyzd1[3];
2213  if(rhit.IsCalibrated()){
2214  gev = rhit.GeV();
2215  }else if(pUseUncalib == true){// use MC calibration if a cell is not calibrated
2216 // std::cout<<"slice, plane, cell, wx, gev "<<ic<<" "<<chit->Plane()<<" "<<chit->Cell()<<" "<<wx<<" Not calibrated?!" << std::endl;
2217  double adc = double(chit->ADC());
2218  double Att = 1;
2219  if(view == geo::kX){
2220  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2221  }else if(view == geo::kY){
2222  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2223  }
2224  gev = adc/Att;
2225  }
2226 
2227  double ecor = 1.;
2228  double disttoelec = geom->Plane(chit->Plane())->Cell(chit->Cell())->HalfL()+pigtail-wx;
2229  if(evt.isRealData()==true&&pCorDataCell)ecor = RecoJMShower::CellEDataMC(gev, disttoelec, view, fDetid); // correct Data/MC difference for cell energy
2230  gev = gev/ecor;
2231  sliceExclGeV += gev;
2232  }
2233  }
2234 
2235  //loop over y cells
2236  for(unsigned int i=0;i<clust->NYCell();i++){
2237  geo::OfflineChan key(clust->YCell(i)->Plane(),clust->YCell(i)->Cell());
2238  //use hits above pe threshold and are not on any showers
2239  if((clust->YCell(i)->PE())>pCellHitPHThresh&&hitshmap[key].size()==0){
2240  //make reco hits out of cells
2241  const art::Ptr<rb::CellHit>& chit = clust->YCell(i);
2242  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, wy ));
2243 
2244  double gev=0;
2245  geo::View_t view = geom->Plane(chit->Plane())->View();
2246  double pigtail = geom->GetPigtail(geom->Plane(chit->Plane())->Cell(chit->Cell())->Id());
2247  //double wxtoelec = wx - (pigtail - 92.6624);//W - (pigtail - avg. pigtail)
2248  //double wytoelec = wy - (pigtail - 92.6624);//W - (pigtail - avg. pigtail)
2249  double w1 = wx;
2250  if (view == geo::kY) w1 = wy;
2251  double xyzd1[4];
2252  calib::GetXYZD(chit->OfflineChan(), w1, xyzd1);
2253  double readDist = xyzd1[3];
2254  if(rhit.IsCalibrated()){
2255  gev = rhit.GeV();
2256  }else if(pUseUncalib == true){// use MC calibration if a cell is not calibrated
2257 // std::cout<<"slice, plane, cell, wy, gev "<<ic<<" "<<chit->Plane()<<" "<<chit->Cell()<<" "<<wy<<" Not calibrated?!" << std::endl;
2258  double adc = double(chit->ADC());
2259  double Att = 1;
2260  if(view == geo::kX){
2261  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2262  }else if(view == geo::kY){
2263  Att = RecoJMShower::CellADCToGeV(readDist, view, fDetid);
2264  }
2265  gev = adc/Att;
2266  }
2267  double ecor = 1.;
2268  double disttoelec = geom->Plane(chit->Plane())->Cell(chit->Cell())->HalfL()+pigtail-wy;
2269  if(evt.isRealData()==true&&pCorDataCell)ecor = RecoJMShower::CellEDataMC(gev, disttoelec, view, fDetid); // correct Data/MC difference for cell energy
2270  gev = gev/ecor;
2271 
2272  sliceExclGeV += gev;
2273  }
2274  }
2275 
2276  //After clustering, initial pclusters
2277  for(unsigned int i = 0; i < vCluster.size(); i++){
2278  std::vector<PCluster> showerpcluster;
2279  std::vector<PCluster> trackpcluster;
2280  showerpcluster.clear();
2281  trackpcluster.clear();
2282  if(vNantag[i]==true){
2283  fShowerPClusterCol.push_back(showerpcluster);
2284  fTrackPClusterCol.push_back(trackpcluster);
2285  continue;
2286  }
2287  if(startplanecol[i]<=stopplanecol[i]){// deal with backward
2288  for(int ip = (int)startplanecol[i];ip<=(int)stopplanecol[i];ip++){
2289  PCluster pcluster;
2290  pcluster.SetPlane(ip);
2291  showerpcluster.push_back(pcluster);
2292  trackpcluster.push_back(pcluster);
2293  }
2294  }else{
2295  for(int ip = (int)startplanecol[i];ip>=(int)stopplanecol[i];ip--){
2296  PCluster pcluster;
2297  pcluster.SetPlane(ip);
2298  showerpcluster.push_back(pcluster);
2299  trackpcluster.push_back(pcluster);
2300  }
2301  }
2302  for (unsigned int nh=0;nh<vCluster[i].NCell();nh++) {
2303  const art::Ptr<rb::CellHit> cellhit = vCluster[i].Cell(nh);
2304  int plane = (int)cellhit->Plane()-(int)startplanecol[i];
2305  if(startplanecol[i]>stopplanecol[i])plane=-plane;
2306  showerpcluster[plane].Add((cellhit));
2307  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
2308  bool trktag=false;// on orignal tracks ?
2309  for(unsigned int imap=0;imap<hittrkmap[key].size();imap++){
2310  if(hittrkmap[key][imap]==(int)i)trktag=true;
2311  }
2312  if(trktag==true)trackpcluster[plane].Add((cellhit));
2313  }
2314  fShowerPClusterCol.push_back(showerpcluster);
2315  fTrackPClusterCol.push_back(trackpcluster);
2316  }
2317  for(unsigned int i = 0; i < vCluster.size(); i++){//calculate vertex doca to each shower
2318  double startx=-9999,starty=-9999,startz=-9999;
2319  double stopx=-9999,stopy=-9999,stopz=-9999;
2320  bool startxtag=false,startytag=false,startztag=false;
2321  bool stopxtag=false,stopytag=false,stopztag=false;
2322  if(vNantag[i]==true){
2323  vtxdocacol.push_back(9999);
2324  continue;
2325  }
2326  for(int plane = 0; plane < (int)fShowerPClusterCol[i].size(); plane++ ){
2327  unsigned int geoplane = fShowerPClusterCol[i][plane].Plane();
2328  geo::View_t view = geom->Plane(geoplane)->View();
2329  unsigned int planecell = RecoJMShower::GetTrkCPlaneCell((unsigned int)plane, i);
2330  if(planecell==9999)continue;
2331  Double_t cellXYZ[3]={0,0,0};
2332  geom->Plane(geoplane)->Cell(planecell)->GetCenter(cellXYZ);
2333  if(startztag==false){startz = cellXYZ[2];startztag=true;}
2334  if(view==geo::kX&&startxtag==false){startx = cellXYZ[0];startxtag=true;}
2335  if(view==geo::kY&&startytag==false){starty = cellXYZ[1];startytag=true;}
2336  if(startxtag==true&&startytag==true&&startztag==true)break;
2337  }
2338 
2339  for(int plane = (int)fShowerPClusterCol[i].size()-1;plane >= 0; plane-- ){
2340  unsigned int geoplane = fShowerPClusterCol[i][plane].Plane();
2341  geo::View_t view = geom->Plane(geoplane)->View();
2342  unsigned int planecell = RecoJMShower::GetTrkCPlaneCell((unsigned int)plane, i);
2343  if(planecell==9999)continue;
2344  Double_t cellXYZ[3]={0,0,0};
2345  geom->Plane(geoplane)->Cell(planecell)->GetCenter(cellXYZ);
2346  if(stopztag==false){stopz = cellXYZ[2];stopztag=true;}
2347  if(view==geo::kX&&stopxtag==false){stopx = cellXYZ[0];stopxtag=true;}
2348  if(view==geo::kY&&stopytag==false){stopy = cellXYZ[1];stopytag=true;}
2349  if(stopxtag==true&&stopytag==true&&stopztag==true)break;
2350  }
2351  double vtxDoca = 9999;
2352  if(!(startx<-9000||starty<-9000||startz<-9000||stopx<-9000||stopy<-9000||stopz<-9000)&&evtvtxcol.size()>0){
2353  TVector3 start(startx,starty,startz);
2354  TVector3 stop(stopx,stopy,stopz);
2355  vtxDoca = RecoJMShower::GetPointDoca(evtvtxcol[0], start, stop);
2356  }
2357  vtxdocacol.push_back(vtxDoca);
2358  }
2359 
2360  //Initial tcclusters
2361  /*
2362  for(unsigned int i = 0; i < vCluster.size(); i++){
2363  std::vector<TCCluster> showertccluster;
2364  showertccluster.clear();
2365  if(vNantag[i]==true){
2366  fShowerTCClusterCol.push_back(showertccluster);
2367  continue;
2368  }
2369  for(int itc = 0; itc<=20;itc++){
2370  TCCluster tccluster;
2371  tccluster.SetTCell(itc);
2372  showertccluster.push_back(tccluster);
2373  }
2374 
2375  for (unsigned int nh=0;nh<vCluster[i].NCell();nh++) {
2376  const art::Ptr<rb::CellHit> cellhit = vCluster[i].Cell(nh);
2377  unsigned transcellId = cellhit->Cell();
2378  unsigned int planecell = RecoJMShower::GetTrkPlaneCell(cellhit->Plane(), vCluster[i]);
2379  int plane = (int)cellhit->Plane()-(int)startplanecol[i];
2380  if(startplanecol[i]>stopplanecol[i])plane=-plane;
2381  // unsigned int planecell = RecoJMShower::GetTrkCPlaneCell((unsigned int)plane, i);
2382  // std::cout<<"GetTrkPlaneCell, GetTrkCPlaneCell "<<planecell1<<" "<<planecell<<std::endl;
2383  if(planecell==9999)continue;
2384  int itc = std::std::abs(int(planecell) - int(transcellId));
2385  if(itc>20)continue;
2386  showertccluster[itc].Add(cellhit);
2387  }
2388  fShowerTCClusterCol.push_back(showertccluster);
2389  }
2390  */
2391 
2392  //fill cells to showercol
2393  for(unsigned int i = 0; i < vCluster.size(); i++){
2394  if(vNantag[i]==true)continue;
2395  // jmshower::JMShower shower(vCluster[i], startcol[i], dircol[i]);
2396  jmshower::JMShower shower(vCluster[i], i, vCluster[i].TotalLength());
2397 
2398  if(vCluster[i].NCell()==0){
2399 
2400  std::cout<<"Zero Cell Track! i, NCell, NTrk "<<i<<" "<<track[i]->NCell()<<" "<<track.size()<<std::endl;
2401  std::cout<<"Zero Cell Shower! i, NCell, NSh "<<i<<" "<<vCluster[i].NCell()<<" "<<vCluster.size()<<std::endl;
2402  continue;
2403  }
2404  // shower.fPClusterCol.clear();
2405  // shower.SetWindowSize(shower.ExtentPlane());
2406  shower.SetWindowSize(fShowerPClusterCol[i].size());
2407  shower.SetStartPlane(startplanecol[i]);
2408  shower.SetStopPlane(stopplanecol[i]);
2409  shower.SetNPlane(std::abs((int)stopplanecol[i]-(int)startplanecol[i])+1);
2410  shower.SetXNPlane(vXCluster[i].ExtentPlane());
2411  shower.SetYNPlane(vYCluster[i].ExtentPlane());
2412  shower.SetXMinPlane(vXCluster[i].MinPlane());
2413  shower.SetYMinPlane(vYCluster[i].MinPlane());
2414  shower.SetXMaxPlane(vXCluster[i].MaxPlane());
2415  shower.SetYMaxPlane(vYCluster[i].MaxPlane());
2416 
2417  shower.SetHitColSize(shower.NCell());
2418  geo::OfflineChan key(vCluster[i].Cell(0)->Plane(),vCluster[i].Cell(0)->Cell());
2419  shower.SetISlice(ic);
2420  shower.SetSliceGeV(sliceGeV);
2421  shower.SetSliceEtot(sliceEtot);
2422  shower.SetSliceEdge10GeV(sliceEdge10GeV);
2423  shower.SetSliceEdge5GeV(sliceEdge5GeV);
2424  shower.SetSliceEdge2GeV(sliceEdge2GeV);
2425 
2426  shower.SetNCellToEdge(shNCellToEdge[i]);
2427  shower.SetClosetEdge(shClosetEdge[i]);
2428  shower.SetEdge20GeV(shEdge20GeV[i]);
2429  shower.SetEdge15GeV(shEdge15GeV[i]);
2430  shower.SetEdge10GeV(shEdge10GeV[i]);
2431  shower.SetEdge5GeV(shEdge5GeV[i]);
2432  shower.SetEdge2GeV(shEdge2GeV[i]);
2433 
2434  shower.SetSliceExclGeV(sliceExclGeV);
2435  for (unsigned int nh=0;nh<shower.NCell();nh++) {
2436  art::Ptr<rb::CellHit> cellhit = shower.Cell(nh);
2437  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
2438  double cellEnergy = fHitEWeightMap[key][i];
2439  shower.SetHitDeconvE(nh, cellEnergy);
2440  shower.SetHitDist(nh, fHitDistMap[key][i]);
2441  shower.SetHitSumE(nh, fHitEMap[key]);
2442  }
2443 
2444  //Set ID
2445  shower.SetID(i);
2446  //Set Start and end point
2447  // shower.AppendTrajectoryPoint(endcol[i]);
2448  shower.SetStop(endcol[i]);
2449  shower.SetStart(startcol[i]);
2450  shower.SetGap(gapcol[i]);
2451 // determine the gap between vertex and the start of the continuous plane
2452  unsigned int newstartp = RecoJMShower::ContStartPlane(shower, 0, 3);
2453  unsigned int contstartp = fShowerPClusterCol[i][newstartp].Plane();
2454  double contGap = 0;
2455  if(evtvtxplanecol.size()>0&&contstartp!=0){
2456  double cellD = 2*geom->Plane(0)->Cell(0)->HalfD();
2457  if(std::abs(shower.Dir().Z())>0.00001)contGap=(double)std::abs((int)contstartp-(int)evtvtxplanecol[0])*cellD/std::abs(shower.Dir().Z());
2458  }
2459  contGap=std::max(gapcol[i],contGap);
2460  shower.SetContGap(contGap);
2461  shower.SetVtxDoca(vtxdocacol[i]);
2462  shower.SetBackTag(backtagcol[i]);
2463 
2464  //Set Geometry
2465  shower.SetIsFiducial(RecoJMShower::IsFiducial(shower));
2466  shower.SetDistToEdge(RecoJMShower::GetTrkDistToEdge(shower));
2467  //Set Energy
2468  double trkcor = RecoJMShower::ECorrMC(trkGeV[i]);
2469  if(trkcor<0.0001||trkcor>1000)trkcor = 1;
2470 
2471  shower.SetDepositEnergy(RecoJMShower::DepositEnergy(shower, shower.ID()));
2472  shower.SetEnergy(RecoJMShower::Energy(shower));// with position correction
2473  sliceEtot+=RecoJMShower::Energy(shower);
2474  // shower.SetEnergy(RecoJMShower::Energy(shower, shower.ID())); // wihout position correction
2475  shower.SetRadius(RecoJMShower::Radius(shower));
2476 
2477 
2478  TVector3 pos((shower.Start()[0]+RecoJMShower::GetShwStop(shower)[0])/2.0,(shower.Start()[1]+RecoJMShower::GetShwStop(shower)[1])/2.0,(shower.Start()[2]+RecoJMShower::GetShwStop(shower)[2])/2.0);
2479  int iReg = 0;
2480  if(pos[0]>=0&&pos[1]>=0)iReg = 0;
2481  if(pos[0]>=0&&pos[1]<0)iReg = 1;
2482  if(pos[0]<0&&pos[1]>=0)iReg = 2;
2483  if(pos[0]<=0&&pos[1]<=0)iReg = 3;
2484 
2485  double gev = RecoJMShower::Energy(shower);
2486 
2487  double elower[11] = {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75};
2488  double eupper[11] = {0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25};
2489  double epoint[11];
2490  for(int ig=0;ig<11;ig++){
2491  epoint[ig] = (elower[ig]+eupper[ig])/2.;
2492  }
2493 
2494  int igev0 = 0;
2495  int igev1 = 0;
2496  if(gev>epoint[10]){igev0=10;igev1=10;}
2497  else if(gev<epoint[0]){igev0=0;igev1=0;}
2498  else{
2499  for(int ig=0;ig<10;ig++){
2500  if(gev>epoint[ig]&&gev<epoint[ig+1]){igev0=ig;igev1=ig+1;}
2501  }
2502  }
2503  if(igev0>10)igev0=10;
2504  if(igev1>10)igev1=10;
2505  double e0 = gev-epoint[igev0];
2506  double e1 = epoint[igev1] - gev;
2507 
2508 
2509  double elowerm[41]= {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75,
2510  5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75,
2511  10.25,10.75,11.25,11.75,12.25,12.75,13.25,13.75,14.25,14.75,
2512  15.25,15.75,16.25,16.75,17.25,17.75,18.25,18.75,19.25,19.75
2513  };
2514 
2515  double eupperm[41]= { 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75,
2516  5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75,
2517  10.25,10.75,11.25,11.75,12.25,12.75,13.25,13.75,14.25,14.75,
2518  15.25,15.75,16.25,16.75,17.25,17.75,18.25,18.75,19.25,19.75, 20.25
2519  };
2520  double epointm[41];
2521  for(int ig=0;ig<41;ig++){
2522  epointm[ig] = (elowerm[ig]+eupperm[ig])/2.;
2523  }
2524 
2525  int igevm0 = 0;
2526  int igevm1 = 0;
2527  if(gev>epointm[40]){igevm0=40;igevm1=40;}
2528  else if(gev<epointm[0]){igevm0=0;igevm1=0;}
2529  else{
2530  for(int ig=0;ig<40;ig++){
2531  if(gev>epointm[ig]&&gev<epointm[ig+1]){igevm0=ig;igevm1=ig+1;}
2532  }
2533  }
2534  if(igevm0>40)igevm0=40;
2535  if(igevm1>40)igevm1=40;
2536  double em0 = gev-epointm[igevm0];
2537  double em1 = epointm[igevm1] - gev;
2538 
2539 
2540  for(unsigned int plane = 0; plane<fShowerPClusterCol[i].size();plane++){
2541  double dedx = RecoJMShower::GetPlaneDedx(shower, plane);
2542  /*
2543  double path = RecoJMShower::GetTrkHitPath(fShowerPClusterCol[i][plane].Plane(), 0, shower);
2544  double ep = RecoJMShower::Energy(fShowerPClusterCol[i][plane], shower.ID());
2545  if(plane==0&&dedx<0.000001){
2546  unsigned int lastplane = fShowerPClusterCol[i].size()-1;
2547  std::cout<<"Ncell, adc, pe "<<fShowerPClusterCol[i][plane].NCell()<<" "<<fShowerPClusterCol[i][plane].TotalADC()<<" "<<fShowerPClusterCol[i][plane].TotalPE()<<" "<<std::endl;
2548  std::cout<<"gev,ep,path "<<fShowerPClusterCol[i][plane].TotalGeV()<<" "<<ep<<" "<<path<<std::endl;
2549  std::cout<<"last Ncell, adc, pe "<<fShowerPClusterCol[i][lastplane].NCell()<<" "<<fShowerPClusterCol[i][lastplane].TotalADC()<<" "<<fShowerPClusterCol[i][lastplane].TotalPE()<<" "<<std::endl;
2550  std::cout<<"lastgev,ep,path "<<fShowerPClusterCol[i][lastplane].TotalGeV()<<" "<<ep<<" "<<path<<std::endl;
2551  }
2552  */
2553 
2554  shower.SetPlaneDedx(plane, dedx);
2555  int planee1cell = RecoJMShower::GetPlaneE1Cell(fShowerPClusterCol[i][plane]);
2556  shower.SetPlaneE1Cell(plane, planee1cell);
2557  int planecentroidcell = RecoJMShower::GetPlaneCentroidCell(fShowerPClusterCol[i][plane],shower);
2558  shower.SetPlaneCentroidCell(plane,planecentroidcell);
2559 
2560  if(pRecoPlaneR==true){
2561 // double planerad = RecoJMShower::Radius(fShowerPClusterCol[i][plane],shower);//Use jmshower
2562  double planerad = RecoJMShower::Radius(fTrackPClusterCol[i][plane],shower);//Use FuzzyKvertex Radius
2563  shower.SetPlaneRadius(plane, planerad);
2564  }
2565  TVector3 hitxyz = RecoJMShower::GetTrkPlanePos(fShowerPClusterCol[i][plane].Plane(), shower);
2566  shower.SetPlaneHitXYZ(plane, hitxyz);
2567 
2568  if(pRecoPID==true){
2569  for( int itype = int(jmshower::kNull)+1; itype != int(jmshower::kElectronSG); ++itype){
2570  double prob = RecoJMShower::GetInterPlaneDedxProb(shower, plane, dedx, itype, pos, iReg, igev0, igev1, e0, e1);
2571  shower.SetPlaneProb(plane, prob, itype);
2572  }
2573 
2574  if(pRecoPID==true){
2575  double prob = RecoJMShower::GetInterPlaneDedxProb(shower, plane, dedx, int(jmshower::kElectronSG), pos, iReg, igevm0, igevm1, em0, em1);
2576  shower.SetPlaneProb(plane, prob, int(jmshower::kElectronSG));
2577  }
2578 
2579 
2580  }
2581  }
2582  for(unsigned int tcell = 0; tcell<20; tcell++){
2583  shower.SetCellTransDedx(tcell, RecoJMShower::GetCellTransDedx(shower, tcell));
2584  }
2585  double sumdr = 0;
2586  double sum10dr = 0;
2587  for(unsigned int plane = 0; plane<fShowerPClusterCol[i].size();plane++){
2588  if((int)plane>=(int)fShowerPClusterCol[i].size()-1-(int)plane)break;
2589  double r1=shower.PlaneRadius(plane);
2590  double r2=shower.PlaneRadius(fShowerPClusterCol[i].size()-1-plane);
2591  double dr=r2-r1;
2592  if(r1>1e-10&&r2>1e-10){
2593  sumdr+=dr;
2594  if(plane<10)sum10dr+=dr;
2595  }
2596  }
2597  shower.SetSumDRadius(sumdr);
2598  shower.SetSum10DRadius(sum10dr);
2599  //Virtual nodes and moments calculation
2600  std::vector<TVector3> nodetocorecol;
2601  std::vector<geo::OfflineChan> nodexidcol;
2602  std::vector<geo::OfflineChan> nodeyidcol;
2603  std::vector<double> nodetocoreenergycol;
2604  nodetocorecol.clear();
2605  nodetocoreenergycol.clear();
2606  nodexidcol.clear();
2607  nodeyidcol.clear();
2608 
2609  if(pRecoNode==true&&shower.TotalLength()<600.){
2610 
2611  double sumE = 0;
2612  TVector3 sumENodeToCore(0,0,0);
2613  double I11 = 0;
2614  double I12 = 0;
2615  double I21 = 0;
2616  double I22 = 0;
2617 
2618  // Using 3-D nodes to calculate shower balance
2619  // TVector3 core(shower.Start()[0],shower.Start()[1],shower.Start()[2]); // original point of shower frame
2620  TVector3 iz = shower.Dir().Unit(); // z in shower frame
2621  TVector3 ix0(1.0,0.0,0);
2622  TVector3 iy0(0.0,1.0,0);
2623  TVector3 ix = iy0.Cross(iz).Unit(); // x axis in shower frame
2624  TVector3 iy = iz.Cross(ix).Unit(); // y axis in shower frame
2625 
2626  for(unsigned int plane = 0; plane<fShowerPClusterCol[i].size()-1;plane++){
2627  if(fShowerPClusterCol[i][plane].NCell()==0)continue;
2628  unsigned int plane1 = fShowerPClusterCol[i][plane].Plane();
2629  if(fShowerPClusterCol[i][plane+1].NCell()==0)continue;
2630  unsigned int plane2 = fShowerPClusterCol[i][plane+1].Plane();
2631  geo::View_t view1 = geom->Plane(plane1)->View();
2632  geo::View_t view2 = geom->Plane(plane2)->View();
2633  if(view1==view2)continue;
2634  TVector3 core = RecoJMShower::GetTrkPlanePos(plane1, shower);
2635  TVector3 pcore = RecoJMShower::GetTrkPlanePos(plane1, shower);
2636  std::vector<TVector3> planenodecol;
2637  std::vector<double> planenodeenergycol;
2638  planenodecol.clear();
2639  planenodeenergycol.clear();
2640  for (unsigned itp1 = 0;itp1<fShowerPClusterCol[i][plane].NCell();itp1++) {
2641  unsigned int cell1 = fShowerPClusterCol[i][plane].Cell(itp1)->Cell();
2642  for(unsigned itp2 = 0;itp2<fShowerPClusterCol[i][plane+1].NCell();itp2++){
2643  unsigned int cell2 = fShowerPClusterCol[i][plane+1].Cell(itp2)->Cell();
2644  TVector3 node = RecoJMShower::GetCellNodePos(plane1, cell1, plane2, cell2);// define 3-D node
2645  double nodez = (pcore - shower.Start()).Mag();
2646  TVector3 nodetopcore = node - pcore;
2647  nodetopcore.SetZ(nodez);
2648  TVector3 nodetocore = node - core;
2649  nodetocore.SetZ(nodez);
2650 
2651 
2652  TVector3 xyz = node - core;//transform xy to shower frame
2653 
2654  // double xp = ix[0]*xyz[0]+iy[0]*xyz[1]+iz[0]*xyz[2];// rotate to shower frame
2655  // double yp = ix[1]*xyz[0]+iy[1]*xyz[1]+iz[1]*xyz[2];
2656  // double zp = ix[2]*xyz[0]+iy[2]*xyz[1]+iz[2]*xyz[2];
2657  // double rxy = sqrt(xp*xp+yp*yp);
2658  // double rxypcore = sqrt(nodetopcore[0]*nodetopcore[0]+nodetopcore[1]*nodetopcore[1]);
2659 
2660  double p = shower.Dir()[0];
2661  double q = shower.Dir()[1];
2662  double r = shower.Dir()[2];
2663  double x0 = xyz[0];
2664  double y0 = xyz[1];
2665  double z0 = xyz[2];
2666  TVector3 vt(q*z0-r*y0,r*x0-p*z0,p*y0-q*x0);
2667  // double angx = vt.Angle(ix);
2668  // double angy = vt.Angle(iy);
2669  // double angz = vt.Angle(iz);
2670 
2671  // double R = sqrt((q*z0-r*y0)*(q*z0-r*y0)+(r*x0-p*z0)*(r*x0-p*z0)+(p*y0-q*x0)*(p*y0-q*x0));// Rotate xy plane
2672  // TVector3 izcross = ix.Cross(vt);
2673  // double angzzcross = izcross.Angle(iz);
2674  // if(std::abs(angzzcross-3.14159265)<0.01)angx = 2*3.14159265-angx;
2675  // double xrp = R*cos(angx);
2676  // double yrp = R*sin(angx);
2677  // double zrp = sqrt(xyz.Mag()*xyz.Mag()-R*R);
2678  /*
2679  std::cout<<"angle z zp "<<angzzcross<<std::endl;
2680  std::cout<<"angx,angy,angz "<<angx<<" "<<angy<<" "<<angz<<std::endl;
2681  std::cout<<"rxy, rxypcore, R "<<rxy<<" "<<rxypcore<<" "<<R<<std::endl;
2682  std::cout<<"Shower Dir "<<shower.Dir()[0]<<" "<<shower.Dir()[1]<<" "<<shower.Dir()[2]<<std::endl;
2683  std::cout<<"ix Dir "<<ix[0]<<" "<<ix[1]<<" "<<ix[2]<<std::endl;
2684  std::cout<<"iy Dir "<<iy[0]<<" "<<iy[1]<<" "<<iy[2]<<std::endl;
2685  std::cout<<"iz Dir "<<iz[0]<<" "<<iz[1]<<" "<<iz[2]<<std::endl;
2686  std::cout<<"xyz in shower Rotat "<<xrp<<" "<<yrp<<" "<<zrp<<std::endl;
2687  std::cout<<"xyz in shower "<<xp<<" "<<yp<<" "<<zp<<std::endl;
2688  std::cout<<"xyz in shower pcore "<<nodetopcore[0]<<" "<<nodetopcore[1]<<" "<<nodez<<std::endl;
2689  */
2690 
2691  // TVector3 nodetocore(xp,yp,zp);
2692  // TVector3 nodetocore(xrp,yrp,zrp);
2693  double rcn = nodetocore.Perp();
2694  geo::OfflineChan key1(plane1,cell1);
2695  geo::OfflineChan key2(plane2,cell2);
2696  double nodeenergy = fHitEWeightMap[key1][i]+fHitEWeightMap[key2][i];
2697  rcn = rcn*nodeenergy;
2698  sumE += nodeenergy;
2699  sumENodeToCore+=nodeenergy*nodetocore;
2700  I11 += nodeenergy*(nodetocore.Y()*nodetocore.Y());
2701  I12 += -nodeenergy*(nodetocore.X()*nodetocore.Y());
2702  I21 += -nodeenergy*(nodetocore.X()*nodetocore.Y());
2703  I22 += nodeenergy*(nodetocore.X()*nodetocore.X());
2704 
2705  planenodecol.push_back(node);
2706  planenodeenergycol.push_back(nodeenergy);
2707 
2708  nodetocorecol.push_back(nodetocore);
2709  nodetocoreenergycol.push_back(nodeenergy);
2710  if(view1 == geo::kX){
2711  nodexidcol.push_back(key1);
2712  nodeyidcol.push_back(key2);
2713  }else{
2714  nodexidcol.push_back(key2);
2715  nodeyidcol.push_back(key1);
2716  }
2717  }
2718  }
2719  }
2720  // calculate eiganvalues and diagonalize inertial tensor
2721  double Ixx = 0.5*(I11+I22+sqrt(pow(I11+I22,2)-4*(I11*I22-I12*I21)));
2722  double Iyy = 0.5*(I11+I22-sqrt(pow(I11+I22,2)-4*(I11*I22-I12*I21)));
2723  double asym = 1;
2724  if(std::abs(Ixx+Iyy)>1e-10)asym=(Ixx-Iyy)/(Ixx+Iyy);
2725  double alpha = 3.14159265/2.;
2726  if(std::abs(I12)>1e-10)alpha = atan((Ixx-I11)/I12);//rotate nodes to principal axes frame
2727  for(int inode=0;inode<(int)nodetocorecol.size();inode++){
2728  nodetocorecol[inode].RotateZ(-alpha);
2729  }
2730 
2731  shower.SetNodes(nodetocorecol, nodetocoreenergycol,nodexidcol,nodeyidcol);
2732  TVector3 ixyz(Ixx,Iyy,0);
2733  shower.SetIxyz(ixyz);
2734  shower.SetAsym(asym);
2735  }else{
2736  shower.SetNodes(nodetocorecol, nodetocoreenergycol,nodexidcol,nodeyidcol);
2737  TVector3 ixyz(0,0,0);
2738  shower.SetIxyz(ixyz);
2739  shower.SetAsym(0);
2740  }
2741 
2742 
2743 
2744  if(pRecoPID==true){
2745  for( int itype = int(jmshower::kNull)+1; itype != int(jmshower::kLast); ++itype)
2746  {
2747  // std::cout<<"in RecoJMShower itype, LLL "<<itype<<" "<<RecoJMShower::GetDedxLongLL(shower, itype)<<std::endl;
2748  // std::cout<<"in RecoJMShower itype, LLT "<<itype<<" "<<RecoJMShower::GetDedxTransLL(shower, itype)<<std::endl;
2749  // std::cout<<"itype Proton = "<<int(jmshower::kProton)<<std::endl;
2750  shower.SetDedxLongLL(RecoJMShower::GetDedxLongLL(shower, itype), itype);
2751  if(pRecoInvLL==true){
2752 // if(itype==int(jmshower::kGamma))shower.SetDedxInvLongLL(RecoJMShower::GetDedxInvLongLL(shower, itype), itype);
2753  shower.SetDedxInvLongLL(RecoJMShower::GetDedxInvLongLL(shower, itype), itype);
2754  }
2755  shower.SetDedxTransLL(RecoJMShower::GetDedxTransLL(shower, itype), itype);
2756  }
2757  }
2758 
2759  shower.SetNMIPPlane(RecoJMShower::NMIPPlane(shower));
2760  shower.SetCentroid(RecoJMShower::GetCentroid(shower));
2761  shower.SetVtxGeV(vtxGeV[i]);
2762  shower.SetVtxCenterGeV(vtxCenterGeV[i]);
2763  shower.SetVtxCenterGeV5cell(vtxCenterGeV5cell[i]);
2764  shower.SetVtxCenterGeV10cell(vtxCenterGeV10cell[i]);
2765  shower.SetVtxCenterGeV15cell(vtxCenterGeV15cell[i]);
2766  shower.SetVtxCenterGeV20cell(vtxCenterGeV20cell[i]);
2767  shower.SetTrkGeV(trkGeV[i]/trkcor);
2768  // for(unsigned int plane=0;plane<fShowerPClusterCol[i].size();plane++){
2769  // shower.fPClusterCol.push_back(fShowerPClusterCol[i][plane]);
2770  // }
2771 
2772  // for(int itc = 0; itc<=20;itc++){
2773  // shower.fTCClusterCol.push_back(fShowerTCClusterCol[i][itc]);
2774  // }
2775 
2776  /*
2777  if(pRecoPID==true)shower.SetIsMuon(RecoJMShower::IsMuon(shower));
2778  if(pRecoANN==1){
2779  shower.SetANN(0, RecoJMShower::CalcANN(0, shower));
2780  std::cout<<"pRecoANN, ANN = "<<pRecoANN<<" "<<RecoJMShower::CalcANN(0, shower)<<std::endl;
2781  }
2782  //std::cout<<"Testing ... "<<_efrac_2slides<<" "<<_efrac_4slides<<" "<<_efrac_6slides<<" "<<_efrac_2sigmaHOR<<" "<<_efrac_4sigmaHOR<<std::endl;
2783  */
2784 
2785 
2786  //fill in the rvp variables that could potentially be used
2787  if(pRecoRVP==true){
2788  RVPInfo rvpstats = this->GetRVPStats(slicecol, ic);
2789  shower.SetEfrac_2slides(rvpstats.efrac_2slides);
2790  shower.SetEfrac_4slides(rvpstats.efrac_4slides);
2791  shower.SetEfrac_6slides(rvpstats.efrac_6slides);
2792  shower.SetEfrac_20slides(rvpstats.efrac_20slides);
2793  shower.SetEfrac_2sigmaHOR(rvpstats.efrac_2sigmaHOR);
2794  shower.SetEfrac_4sigmaHOR(rvpstats.efrac_4sigmaHOR);
2795  shower.SetEiso_2sigmaHOR(rvpstats.eiso_2sigmaHOR);
2796  shower.SetEiso_4sigmaHOR(rvpstats.eiso_4sigmaHOR);
2797  shower.SetEiso_6sigmaHOR(rvpstats.eiso_6sigmaHOR);
2798  }
2799  showercol.push_back(shower);
2800  }
2801 // std::cout<<"n sh "<<showercol.size()<<std::endl;
2802  std::vector<TLorentzVector> showerpGam;// Set shower momentum
2803  showerpGam.clear();
2804 
2805  for(unsigned int i = 0; i < showercol.size(); i++){
2806  double showereraw = showercol[i].Energy();
2807  TLorentzVector showerptrk(showercol[i].Dir()[0]*showereraw,showercol[i].Dir()[1]*showereraw,showercol[i].Dir()[2]*showereraw,showereraw);
2808  showerpGam.push_back(showerptrk);
2809  }
2810 
2811  for(unsigned int i = 0; i < showercol.size(); i++){
2812  double minpi0mgg = -9999;
2813  int minj = -1;
2814  double shexcle = 0;
2815  for(unsigned int j=0;j< showercol.size();j++){ // pi0 mass calculation
2816  if(j==i)continue;
2817  if(showercol[j].ISlice()!=showercol[i].ISlice())continue;
2818  TLorentzVector p2g = showerpGam[i] + showerpGam[j];
2819  double mgg = p2g.M();
2820  if( std::abs(minpi0mgg-0.15) > std::abs(mgg-0.15)){
2821  minj = (int)j;
2822  minpi0mgg = mgg;
2823  }
2824  shexcle+=showercol[j].DepositEnergy();
2825  }
2826  shexcle+=sliceExclGeV;
2827  if(minpi0mgg<0)minpi0mgg=0.;
2828  showercol[i].SetPi0Mgg(minpi0mgg);
2829  showercol[i].SetPi0G2ID(minj);
2830  showercol[i].SetSliceEtot(sliceEtot);
2831  showercol[i].SetExclEnergy(shexcle);
2832  double elecE = showercol[i].Energy();
2833  double hadE = showercol[i].ExclEnergy();
2834  //Nue Energy = electron energy + hadron energy
2835  double nuE = (elecE+hadE/(4.00599e-01
2836  -5.08688e+00*TMath::Exp(-1.07930e+01*hadE-(-1.17127e+01)*pow(hadE,2)-1.01452e+01*pow(hadE,3)) +4.71995e+00*TMath::Exp(-1.02536e+01*hadE-(-9.81034e+00)*pow(hadE,2)-4.22648e+00*pow(hadE,3))));
2837  // double nuE = (showercol[i].Energy() + 0.282525 + 1.0766*(hadE))/1.015;
2838  showercol[i].SetNueEnergy(nuE);
2839  }
2840 
2841 
2842  // TVector3 nuDir0(-0.845e-3, 52.552e-3, 0.998618);
2843  TVector3 nuDir = geom->NuMIBeamDirection();
2844  for(unsigned int i = 0; i < showercol.size(); i++){
2845  double Theta = showerpGam[i].Angle(nuDir);
2846  showercol[i].SetTheta(Theta);
2847  }
2848 
2849  for(unsigned int i=0;i < showercol.size(); i++){
2850  if(pRecoPID==true){
2851  showercol[i].SetIsMuon(RecoJMShower::IsMuon(showercol[i]));
2852  showercol[i].SetIsInvPhoton(RecoJMShower::IsInvPhoton(showercol[i]));
2853  }
2854  if(pRecoANN==1){
2855  showercol[i].SetANN(0, RecoJMShower::CalcANN(0, showercol[i]));
2856  showercol[i].SetANN(2, RecoJMShower::CalcANN(2, showercol[i]));
2857  showercol[i].SetANN(3, RecoJMShower::CalcANN(3, showercol[i]));
2858  showercol[i].SetANN(4, RecoJMShower::CalcANN(4, showercol[i]));
2859  showercol[i].SetANN(5, RecoJMShower::CalcANN(5, showercol[i]));
2860  showercol[i].SetANN(6, RecoJMShower::CalcANN(6, showercol[i]));
2861  showercol[i].SetANN(7, RecoJMShower::CalcANN(7, showercol[i]));
2862  // showercol[i].SetANN(1, RecoJMShower::CalcANN(1, showercol[i]));
2863  // std::cout<<"pRecoANN, ANN = "<<pRecoANN<<" "<<RecoJMShower::CalcANN(0, showercol[i])<<std::endl;
2864  }
2865  }
2866 
2867  for(unsigned int i = 0;i<showercol.size();i++){
2868  recoshowercol->push_back(showercol[i]);
2869  evtshowercol.push_back(showercol[i]);
2870 //ad a protection for reconstruction fail
2871  if(i>=fprongidx.size())continue;
2872  if((unsigned int)fprongidx[i]>=track.size())continue;
2873  util::CreateAssn(prod,evt, *recoshowercol, track[fprongidx[i]] ,*prongassns);
2874  }
2875  }
2876 
2877 
2878 
2879  return evtshowercol;
2880  }
2881 
2882  std::vector<jmshower::JMShower> RecoJMShower::ExclShowers(const art::Event& evt){
2883  return fExcludeShowerCol;
2884  }
2885 
2886  bool RecoJMShower::PreSel(const art::Event& evt){
2887  bool tagPreSel = true;
2888  art::Handle< std::vector<rb::Prong> > trackhandle;
2889  std::vector<art::Ptr<rb::Prong> > track;
2890  evt.getByLabel(pTrackDir, trackhandle);
2891 
2892  for(unsigned int i = 0; i < trackhandle->size(); i++){
2893  art::Ptr<rb::Prong> p(trackhandle, i);
2894  track.push_back(p);
2895  }
2896 
2897  for(unsigned int itrk=0;itrk<track.size();itrk++){
2898  if(track[itrk]->ExtentPlane()>200)tagPreSel = false;
2899  }
2900  return tagPreSel;
2901  }
2902 
2903 
2904 //******************************************************
2905 // Geoometry
2906 //******************************************************
2907 
2908  double RecoJMShower::GetPointDoca(TVector3 vtx, TVector3 start, TVector3 stop){
2909  double d = 9999;
2910  double denorm = (start-stop).Mag();
2911  if(denorm<0.00001){d = (vtx-start).Mag();}else{
2912  d = ((vtx-start).Cross(vtx-stop)).Mag()/denorm;
2913  TVector3 dir = stop-start;
2914  }
2915  return d;
2916  }
2917 
2918 
2919  TVector3 RecoJMShower::GetCellDistToPoint(unsigned int plane, unsigned int cell, TVector3 point){
2920  Double_t cellXYZ[3]={0,0,0};
2921  art::ServiceHandle<geo::Geometry> geom;
2922  geom->Plane(plane)->Cell(cell)->GetCenter(cellXYZ);
2923  geo::View_t view = geom->Plane(plane)->View();
2924  double celldistx = -99;
2925  double celldisty = -99;
2926  double celldistz = -99;
2927  if(view == geo::kX){
2928  celldistx = std::abs(cellXYZ[0]-point[0]);
2929  celldisty = 0;
2930  celldistz = std::abs(cellXYZ[2]-point[2]);
2931  }else if(view == geo::kY){
2932  celldistx = 0;
2933  celldisty = std::abs(cellXYZ[1]-point[1]);
2934  celldistz = std::abs(cellXYZ[2]-point[2]);
2935  }
2936  TVector3 celldist(celldistx,celldisty,celldistz);
2937  return celldist;
2938  }
2939 
2940  //******************************************************************************************************************
2941  double RecoJMShower::GetCellDistToTrk(unsigned int plane, unsigned int cell, TVector3 start, TVector3 stop){
2942  Double_t cellXYZ[3]={0,0,0};
2943  art::ServiceHandle<geo::Geometry> geom;
2944  geom->Plane(plane)->Cell(cell)->GetCenter(cellXYZ);
2945  geo::View_t view = geom->Plane(plane)->View();
2946 
2947  double x1 = start[0];
2948  double y1 = start[1];
2949  double z1 = start[2];
2950  double p1 = (stop[0] - start[0])/(stop - start).Mag();
2951  double q1 = (stop[1] - start[1])/(stop - start).Mag();
2952  double r1 = (stop[2] - start[2])/(stop - start).Mag();
2953 
2954  double x2 = cellXYZ[0];
2955  double y2 = cellXYZ[1];
2956  double z2 = cellXYZ[2];
2957 
2958  double p2 = 1;
2959  double q2 = 0;
2960  double r2 = 0;
2961  if(view == geo::kX){
2962  p2 = 0;
2963  q2 = 1;
2964  r2 = 0;
2965  }else if(view == geo::kY){
2966  p2 = 1;
2967  q2 = 0;
2968  r2 = 0;
2969  }
2970  double ma = x2-x1;
2971  double mb = y2-y1;
2972  double mc = z2-z1;
2973  double md = p1;
2974  double me = q1;
2975  double mf = r1;
2976  double mg = p2;
2977  double mh = q2;
2978  double mi = r2;
2979 
2980  TMatrixD dDeNorm1(2,2);
2981  TMatrixD dDeNorm2(2,2);
2982  TMatrixD dDeNorm3(2,2);
2983  dDeNorm1[0][0] = p1;
2984  dDeNorm1[0][1] = q1;
2985  dDeNorm1[1][0] = p2;
2986  dDeNorm1[1][1] = q2;
2987  double det1 = dDeNorm1[0][0]*dDeNorm1[1][1] - dDeNorm1[0][1]*dDeNorm1[1][0];
2988  dDeNorm2[0][0] = q1;
2989  dDeNorm2[0][1] = r1;
2990  dDeNorm2[1][0] = q2;
2991  dDeNorm2[1][1] = r2;
2992  double det2 = dDeNorm2[0][0]*dDeNorm2[1][1] - dDeNorm2[0][1]*dDeNorm2[1][0];
2993  dDeNorm3[0][0] = r1;
2994  dDeNorm3[0][1] = p1;
2995  dDeNorm3[1][0] = r2;
2996  dDeNorm3[1][1] = p2;
2997  double det3 = dDeNorm3[0][0]*dDeNorm3[1][1] - dDeNorm3[0][1]*dDeNorm3[1][0];
2998  double det0 = ma*me*mi + mb*mf*mg + mc*md*mh - mc*me*mg - mb*md*mi - ma*mf*mh;
2999  double celldist = std::abs(det0)/sqrt(det1*det1+det2*det2+det3*det3);
3000  return celldist;
3001  }
3002 
3003 //****************************************************************************************************************
3004 
3005 
3006  TVector3 RecoJMShower::GetShwStop(rb::Prong shw){
3007  TVector3 shwdir=shw.Dir();
3008  TVector3 shwstart=shw.Start();
3009  TVector3 shwend = shwstart + shw.TotalLength()*shwdir;
3010  return shwend;
3011  }
3012 
3013  TVector3 RecoJMShower::GetShwStop(art::Ptr<rb::Prong> shw){
3014  TVector3 shwdir=shw->Dir();
3015  TVector3 shwstart=shw->Start();
3016  TVector3 shwend = shwstart + shw->TotalLength()*shwdir;
3017  return shwend;
3018  }
3019 
3020 
3021 //given a plane, cell and track, return distance of cell with respect to track
3022  TVector3 RecoJMShower::GetTrkHitPos(unsigned int plane, unsigned int cell, art::Ptr<rb::Prong> trk){
3023  Double_t cellXYZ[3]={0,0,0}; //declare container to hold cell coordinate
3024  art::ServiceHandle<geo::Geometry> geom;
3025  geom->Plane(plane)->Cell(cell)->GetCenter(cellXYZ); //get cell center
3026 
3027  //get track direction, start, end
3028  TVector3 trkdir=trk->Dir();
3029  TVector3 trkstart=trk->Start();
3030  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3031 // double cellL=2*geom->Plane(plane)->Cell(cell)->HalfL();
3032  geo::View_t view = geom->Plane(plane)->View();
3033  double celldist = 0;
3034  double trkx = 0;
3035  double trky = 0;
3036  //calculate the tracks x position at the cell's z coordinate
3037  trkx = trkstart.x()+(trkdir.x()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3038  //calculate the tracks y position at the cell's z coordinate
3039  trky = trkstart.y()+(trkdir.y()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3040 // double trkdtr0 = 0;
3041  double trkdtr = 0;
3042  if(view == geo::kX){
3043  //calculate distance from cell x coordinate to track x coordninate at same z position
3044  //NOTE: Perpendicular distance to track may be a better metric as this will be incorrect for steep tracks
3045  celldist=std::abs((cellXYZ[0]-trkstart.x())-((cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z())));
3046 // trkdtr0 = geom->Plane(plane)->Cell(cell)->DistToReadOut(trky); //distance to readout for cell
3047  trkdtr = geom->DistToElectronics(trky, *geom->Plane(plane)->Cell(cell));
3048 
3049  }else if(view == geo::kY){
3050  //calculate distance from cell x coordinate to track x coordninate at same z position
3051  celldist=std::abs((cellXYZ[1]-trkstart.y())-((cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z())));
3052 // trkdtr0 = geom->Plane(plane)->Cell(cell)->DistToReadOut(trkx); //distance to readout for cell
3053  trkdtr = geom->DistToElectronics(trkx, *geom->Plane(plane)->Cell(cell));
3054  }
3055  if(trkdtr<0)trkdtr=0;
3056  //return results
3057  TVector3 def(-1000,-1000,-1000);
3058  TVector3 trkhitposxcell(trkdtr,trky,celldist);
3059  TVector3 trkhitposycell(trkdtr,trkx,celldist);
3060  if(view == geo::kX){return trkhitposxcell;}
3061  if(view == geo::kY){return trkhitposycell;}
3062  return def;
3063  }
3064 
3065 //*********************************************************************************************************************
3066  TVector3 RecoJMShower::GetTrkHitPos(unsigned int plane, unsigned int cell, rb::Prong trk){
3067  Double_t cellXYZ[3]={0,0,0};
3068  art::ServiceHandle<geo::Geometry> geom;
3069  geom->Plane(plane)->Cell(cell)->GetCenter(cellXYZ);
3070  TVector3 trkdir=trk.Dir();
3071  TVector3 trkstart=trk.Start();
3072  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3073 // double cellL=2*geom->Plane(plane)->Cell(cell)->HalfL();
3074  geo::View_t view = geom->Plane(plane)->View();
3075  double celldist = 0;
3076  double trkx = 0;
3077  double trky = 0;
3078  trkx = trkstart.x()+(trkdir.x()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3079  trky = trkstart.y()+(trkdir.y()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3080 
3081 // double trkdtr0 = 0;
3082  double trkdtr = 0;
3083  if(view == geo::kX){
3084  celldist=std::abs((cellXYZ[0]-trkstart.x())-((cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z())));
3085 // trkdtr = geom->Plane(plane)->Cell(cell)->DistToReadOut(trky);
3086  trkdtr = geom->DistToElectronics(trky, *geom->Plane(plane)->Cell(cell));
3087  }else if(view == geo::kY){
3088  celldist=std::abs((cellXYZ[1]-trkstart.y())-((cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z())));
3089 // trkdtr = geom->Plane(plane)->Cell(cell)->DistToReadOut(trkx);
3090  trkdtr = geom->DistToElectronics(trkx, *geom->Plane(plane)->Cell(cell));
3091  }
3092  if(trkdtr<0)trkdtr=0;
3093  TVector3 def(-1000,-1000,-1000);
3094  TVector3 trkhitposxcell(trkdtr,trky,celldist);
3095  TVector3 trkhitposycell(trkdtr,trkx,celldist);
3096  if(view == geo::kX){return trkhitposxcell;}
3097  if(view == geo::kY){return trkhitposycell;}
3098  return def;
3099  }
3100 
3101  TVector3 RecoJMShower::GetTrkHitPos(unsigned int plane, unsigned int cell, TVector3 trkdir, TVector3 trkstart, TVector3 trkstop){
3102  Double_t cellXYZ[3]={0,0,0}; // XYZ[0] = x, XYZ[1] = y, XYZ[2] = z
3103  art::ServiceHandle<geo::Geometry> geom;
3104  geom->Plane(plane)->Cell(cell)->GetCenter(cellXYZ);
3105 // double cellL=2*geom->Plane(plane)->Cell(cell)->HalfL();
3106  geo::View_t view = geom->Plane(plane)->View();
3107  double celldist = 0;
3108  double trkx = 0;
3109  double trky = 0;
3110  trkx = trkstart.x()+(trkdir.x()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3111  trky = trkstart.y()+(trkdir.y()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3112 // double trkdtr0 = 0;
3113  double trkdtr = 0;
3114  if(view == geo::kX){
3115  celldist=std::abs((cellXYZ[0]-trkstart.x())-((cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z())));
3116 // trkdtr0 = geom->Plane(plane)->Cell(cell)->DistToReadOut(trky);
3117  trkdtr = geom->DistToElectronics(trky, *geom->Plane(plane)->Cell(cell));
3118  }else if(view == geo::kY){
3119  celldist=std::abs((cellXYZ[1]-trkstart.y())-((cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z())));
3120 // trkdtr0 = geom->Plane(plane)->Cell(cell)->DistToReadOut(trkx);
3121  trkdtr = geom->DistToElectronics(trkx, *geom->Plane(plane)->Cell(cell));
3122  }
3123  if(trkdtr<0)trkdtr=0;
3124  TVector3 def(-1000,-1000,-1000);
3125  TVector3 trkhitposxcell(trkdtr,trky,celldist);
3126  TVector3 trkhitposycell(trkdtr,trkx,celldist);
3127  if(view == geo::kX){return trkhitposxcell;}
3128  if(view == geo::kY){return trkhitposycell;}
3129  return def;
3130  }
3131 
3132 
3133  TVector3 RecoJMShower::GetTrkPlanePos(unsigned int plane, rb::Prong trk){
3134  Double_t cellXYZ[3]={0,0,0};
3135  art::ServiceHandle<geo::Geometry> geom;
3136  geom->Plane(plane)->Cell(0)->GetCenter(cellXYZ);
3137  TVector3 trkdir=trk.Dir();
3138  TVector3 trkstart=trk.Start();
3139  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3140  double trkplanex=99999,trkplaney=99999, trkplanez=0;
3141  if(std::abs(trkdir.z())>0.0001){
3142  trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3143  trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3144  }else{
3145  trkplanex = (trkstart.x()+trkstop.x())/2.0;
3146  trkplaney = (trkstart.y()+trkstop.y())/2.0;
3147  }
3148  trkplanez = cellXYZ[2];
3149  TVector3 trkplanepos(trkplanex, trkplaney, trkplanez);
3150  return trkplanepos;;
3151  }
3152 
3153 /*
3154  unsigned int RecoJMShower::GetTrkPlaneCell(unsigned int plane, art::Ptr<rb::Prong> trk){
3155  Double_t cellXYZ[3]={0,0,0};
3156  art::ServiceHandle<geo::Geometry> geom;
3157  geom->Plane(plane)->Cell(0)->GetCenter(cellXYZ);
3158  double cellW = 2*geom->Plane(plane)->Cell(0)->HalfW();
3159  TVector3 trkdir=trk->Dir();
3160  TVector3 trkstart=trk->Start();
3161  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3162  double trkplanex=99999,trkplaney=99999, trkplanez=0;
3163  if(std::abs(trkdir.z())>0.0001){
3164  trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3165  trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3166  }
3167  trkplanez = cellXYZ[2];
3168  geo::View_t view = geom->Plane(plane)->View();
3169  int trkcellId = -9999;
3170  if(view == geo::kX)trkcellId = int((trkplanex+cellW/2-cellXYZ[0])/cellW);
3171  if(view == geo::kY)trkcellId = int((trkplaney+cellW/2-cellXYZ[1])/cellW);
3172  return unsigned(trkcellId);
3173  }
3174 */
3175 
3176 
3177  unsigned int RecoJMShower::GetTrkPlaneCell(unsigned int plane, art::Ptr<rb::Prong> trk){
3178 
3179  Double_t cellXYZ[3]={0,0,0};
3180  Double_t cellXYZ1[3];
3181  Double_t cellXYZmax[3];
3182  art::ServiceHandle<geo::Geometry> geom;
3183  geom->Plane(plane)->Cell(0)->GetCenter(cellXYZ);
3184  geom->Plane(plane)->Cell(geom->Plane(plane)->Ncells()-1)->GetCenter(cellXYZ1);
3185  geom->Plane(plane)->Cell(geom->Plane(plane)->Ncells()-1)->GetCenter(cellXYZmax);
3186  double cellW = 2*geom->Plane(plane)->Cell(0)->HalfW();
3187  double cellwx = (cellXYZ1[0] - cellXYZ[0] + cellW)/geom->Plane(plane)->Ncells();
3188  double cellwy = (cellXYZ1[1] - cellXYZ[1] + cellW)/geom->Plane(plane)->Ncells();
3189  TVector3 trkdir=trk->Dir();
3190  TVector3 trkstart=trk->Start();
3191  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3192  double trkplanex=99999,trkplaney=99999;//, trkplanez=0;
3193 // double trkplanex0=99999,trkplaney0=99999;
3194 // double trkplanex1=99999,trkplaney1=99999;
3195  if(std::abs(trkdir.z())>0.0001){
3196  trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3197  trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3198 // trkplanex1 = cellXYZmax[0];
3199 // trkplaney1 = cellXYZmax[1];
3200  }else{
3201  trkplanex = (trkstart.x()+ trkstop.x())/2.0;
3202  trkplaney = (trkstart.y()+ trkstop.y())/2.0;
3203  }
3204 // trkplanez = cellXYZ[2];
3205  geo::View_t view = geom->Plane(plane)->View();
3206 /*
3207  if(view == geo::kX){
3208  std::cout<<"x view plane"<<std::endl;
3209  std::cout<<"x1-x0,cellW "<<" "<<cellXYZ1[0] - cellXYZ[0]<<" "<<cellW<<std::endl;
3210  }else if(view == geo::kY){
3211  std::cout<<"y view plane"<<std::endl;
3212  std::cout<<"y1-y0,cellW "<<" "<<cellXYZ1[1] - cellXYZ[1]<<" "<<cellW<<std::endl;
3213  }
3214 */
3215  int trkcellId = -9999;
3216 // int trkcellId0 = -9999;
3217 // int trkcellId1 = -9999;
3218 // int trkcellIdmax = geom->Plane(plane)->Ncells()-1;
3219  if(view == geo::kX)trkcellId = int((trkplanex+cellwx/2-cellXYZ[0])/cellwx);
3220  if(view == geo::kY)trkcellId = int((trkplaney+cellwy/2-cellXYZ[1])/cellwy);
3221 // if(view == geo::kX)trkcellId0 = int((trkplanex0+cellwx/2-cellXYZ[0])/cellwx);
3222 // if(view == geo::kY)trkcellId0 = int((trkplaney0+cellwy/2-cellXYZ[1])/cellwy);
3223 // if(view == geo::kX)trkcellId1 = int((trkplanex1+cellwx/2-cellXYZ[0])/cellwx);
3224 // if(view == geo::kY)trkcellId1 = int((trkplaney1+cellwy/2-cellXYZ[1])/cellwy);
3225  return unsigned(trkcellId);
3226  }
3227 
3228  unsigned int RecoJMShower::GetTrkPlaneCell(unsigned int plane, rb::Prong trk){
3229 
3230  Double_t cellXYZ[3]={0,0,0};
3231  Double_t cellXYZ1[3];
3232  Double_t cellXYZmax[3];
3233  art::ServiceHandle<geo::Geometry> geom;
3234  geom->Plane(plane)->Cell(0)->GetCenter(cellXYZ);
3235  geom->Plane(plane)->Cell(geom->Plane(plane)->Ncells()-1)->GetCenter(cellXYZ1);
3236  geom->Plane(plane)->Cell(geom->Plane(plane)->Ncells()-1)->GetCenter(cellXYZmax);
3237  double cellW = 2*geom->Plane(plane)->Cell(0)->HalfW();
3238  double cellwx = (cellXYZ1[0] - cellXYZ[0] + cellW)/geom->Plane(plane)->Ncells();
3239  double cellwy = (cellXYZ1[1] - cellXYZ[1] + cellW)/geom->Plane(plane)->Ncells();
3240  TVector3 trkdir=trk.Dir();
3241  TVector3 trkstart=trk.Start();
3242  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3243  double trkplanex=99999,trkplaney=99999;//, trkplanez=0;
3244 // double trkplanex0=99999,trkplaney0=99999;
3245 // double trkplanex1=99999,trkplaney1=99999;
3246  if(std::abs(trkdir.z())>0.0001){
3247  trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3248  trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3249 // trkplanex0 = cellXYZ[0];
3250 // trkplaney0 = cellXYZ[1];
3251 // trkplanex1 = cellXYZmax[0];
3252 // trkplaney1 = cellXYZmax[1];
3253  }else{
3254  trkplanex = (trkstart.x()+trkstop.x())/2.0;
3255  trkplaney = (trkstart.y()+trkstop.y())/2.0;
3256  }
3257 // trkplanez = cellXYZ[2];
3258  geo::View_t view = geom->Plane(plane)->View();
3259 
3260 /*
3261  if(view == geo::kX){
3262  std::cout<<"x view plane"<<std::endl;
3263  std::cout<<"x1-x0,cellW "<<" "<<cellXYZ1[0] - cellXYZ[0]<<" "<<cellW<<std::endl;
3264  }else if(view == geo::kY){
3265  std::cout<<"y view plane"<<std::endl;
3266  std::cout<<"y1-y0,cellW "<<" "<<cellXYZ1[1] - cellXYZ[1]<<" "<<cellW<<std::endl;
3267  }
3268 */
3269 
3270  int trkcellId = -9999;
3271 // int trkcellId0 = -9999;
3272 // int trkcellId1 = -9999;
3273 // int trkcellIdmax = geom->Plane(plane)->Ncells()-1;
3274 
3275  if(view == geo::kX)trkcellId = int((trkplanex+cellwx/2-cellXYZ[0])/cellwx);
3276  if(view == geo::kY)trkcellId = int((trkplaney+cellwy/2-cellXYZ[1])/cellwy);
3277 
3278 // if(view == geo::kX)trkcellId0 = int((trkplanex0+cellwx/2-cellXYZ[0])/cellwx);
3279 // if(view == geo::kY)trkcellId0 = int((trkplaney0+cellwy/2-cellXYZ[1])/cellwy);
3280 //
3281 // if(view == geo::kX)trkcellId1 = int((trkplanex1+cellwx/2-cellXYZ[0])/cellwx);
3282 // if(view == geo::kY)trkcellId1 = int((trkplaney1+cellwy/2-cellXYZ[1])/cellwy);
3283  return unsigned(trkcellId);
3284  }
3285 
3286  unsigned int RecoJMShower::GetTrkCPlaneCell(unsigned int plane, unsigned int i){//calculate centroid of a plane
3287  art::ServiceHandle<geo::Geometry> geom;
3288  double sumid = 0;
3289  double gev = 0;
3290  for (unsigned int nh=0;nh<fShowerPClusterCol[i][plane].NCell();nh++) {
3291  art::Ptr<rb::CellHit> cellhit = fShowerPClusterCol[i][plane].Cell(nh);
3292  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
3293  double cellEnergy = fHitEWeightMap[key][i];
3294  double cellPos = (double)cellhit->Cell() + 0.5;
3295  sumid += cellEnergy*cellPos;
3296  gev += cellEnergy;
3297  }
3298  if(gev<0.00001)return 9999;
3299  double avgid = sumid/gev;
3300  int cid = (int)avgid;
3301  if(cid<0)cid=0;
3302  if(cid>(int)geom->Plane(plane)->Ncells()-1)cid=(int)geom->Plane(plane)->Ncells()-1;
3303  return (unsigned int)cid;
3304  }
3305 
3306  TVector3 RecoJMShower::GetTrkPlaneDistToEdge(unsigned int plane, art::Ptr<rb::Prong> trk){
3307  Double_t cellXYZ[3]={0,0,0};
3308  art::ServiceHandle<geo::Geometry> geom;
3309  geom->Plane(plane)->Cell(0)->GetCenter(cellXYZ);
3310  TVector3 trkdir=trk->Dir();
3311  TVector3 trkstart=trk->Start();
3312  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3313  double trkplanex=99999,trkplaney=99999, trkplanez=0;
3314  if(std::abs(trkdir.z())>0.0001){
3315  trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3316  trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3317  }
3318  trkplanez = cellXYZ[2];
3319 
3320  double xe[4]={9999,9999,9999,9999};
3321  double ye[4]={9999,9999,9999,9999};
3322  double ze[4]={9999,9999,9999,9999};
3323 
3324  xe[0] = std::abs(pXEdge1-trkplanex);
3325  xe[1] = std::abs(pXEdge2-trkplanex);
3326  xe[2] = std::abs(pXEdge3-trkplanex);
3327  xe[3] = std::abs(pXEdge4-trkplanex);
3328 
3329  ye[0] = std::abs(pYEdge1-trkplaney);
3330  ye[1] = std::abs(pYEdge2-trkplaney);
3331  ye[2] = std::abs(pYEdge3-trkplaney);
3332  ye[3] = std::abs(pYEdge4-trkplaney);
3333 
3334  ze[0] = std::abs(pZEdge1-trkplanez);
3335  ze[1] = std::abs(pZEdge2-trkplanez);
3336  ze[2] = std::abs(pZEdge3-trkplanez);
3337  ze[3] = std::abs(pZEdge4-trkplanez);
3338 
3339  double xtoedge = 9999.;
3340  double ytoedge = 9999.;
3341  double ztoedge = 9999.;
3342 
3343  for(int ii=0;ii<4;ii++){
3344  if(xtoedge>xe[ii])xtoedge=xe[ii];
3345  if(ytoedge>ye[ii])ytoedge=ye[ii];
3346  if(ztoedge>ze[ii])ztoedge=ze[ii];
3347  }
3348 
3349 // double xtoedge0 = std::abs(geom->DetHalfWidth()-std::abs(trkplanex));
3350 // double ytoedge0 = std::abs(geom->DetHalfHeight()-std::abs(trkplaney));
3351 // double ztoedge0 = std::min(std::abs(geom->DetLength()-std::abs(trkplanez)),trkplanez);
3352 // std::cout<<"Det W/2,H/2,L "<<geom->DetHalfWidth()<<" "<<geom->DetHalfHeight()<<" "<<geom->DetLength()<<std::endl;
3353 // double cellL = 2*geom->Plane(plane)->Cell(0)->HalfL();
3354 // double cellW = 2*geom->Plane(plane)->Cell(0)->HalfW();
3355 // double cellD = 2*geom->Plane(plane)->Cell(0)->HalfD();
3356 // geo::View_t view = geom->Plane(plane)->View();
3357 /*
3358  if(view == geo::kX){
3359  std::cout<<"x view plane"<<std::endl;
3360  }else if(view == geo::kY){
3361  std::cout<<"y view plane"<<std::endl;
3362  }
3363  std::cout<<"cell L, cellW, cell D "<<cellL<<" "<<cellW<<" "<<cellD<<std::endl;
3364 */
3365 // std::cout<<"To Edge X,Y,Z "<<xtoedge<<" "<<ytoedge<<" "<<ztoedge<<std::endl;
3366  TVector3 toedge(xtoedge, ytoedge, ztoedge);
3367  return toedge;
3368  }
3369 
3370 
3371  TVector3 RecoJMShower::GetTrkPlaneDistToEdge(unsigned int plane, rb::Prong trk){
3372  Double_t cellXYZ[3]={0,0,0};
3373  art::ServiceHandle<geo::Geometry> geom;
3374  geom->Plane(plane)->Cell(0)->GetCenter(cellXYZ);
3375  TVector3 trkdir=trk.Dir();
3376  TVector3 trkstart=trk.Start();
3377  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3378  double trkplanex=99999,trkplaney=99999, trkplanez=0;
3379  if(std::abs(trkdir.z())>0.0001){
3380  trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3381  trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3382  }
3383  trkplanez = cellXYZ[2];
3384 
3385  double xe[4]={9999,9999,9999,9999};
3386  double ye[4]={9999,9999,9999,9999};
3387  double ze[4]={9999,9999,9999,9999};
3388 
3389  xe[0] = std::abs(pXEdge1-trkplanex);
3390  xe[1] = std::abs(pXEdge2-trkplanex);
3391  xe[2] = std::abs(pXEdge3-trkplanex);
3392  xe[3] = std::abs(pXEdge4-trkplanex);
3393 
3394  ye[0] = std::abs(pYEdge1-trkplaney);
3395  ye[1] = std::abs(pYEdge2-trkplaney);
3396  ye[2] = std::abs(pYEdge3-trkplaney);
3397  ye[3] = std::abs(pYEdge4-trkplaney);
3398 
3399  ze[0] = std::abs(pZEdge1-trkplanez);
3400  ze[1] = std::abs(pZEdge2-trkplanez);
3401  ze[2] = std::abs(pZEdge3-trkplanez);
3402  ze[3] = std::abs(pZEdge4-trkplanez);
3403 
3404  double xtoedge = 9999.;
3405  double ytoedge = 9999.;
3406  double ztoedge = 9999.;
3407 
3408  for(int ii=0;ii<4;ii++){
3409  if(xtoedge>xe[ii])xtoedge=xe[ii];
3410  if(ytoedge>ye[ii])ytoedge=ye[ii];
3411  if(ztoedge>ze[ii])ztoedge=ze[ii];
3412  }
3413 
3414 // double xtoedge0 = std::abs(geom->DetHalfWidth()-std::abs(trkplanex));
3415 // double ytoedge0 = std::abs(geom->DetHalfHeight()-std::abs(trkplaney));
3416 // double ztoedge0 = std::min(std::abs(geom->DetLength()-std::abs(trkplanez)),trkplanez);
3417 
3418 // double xtoedge = geom->DetHalfWidth()-std::abs(trkplanex);
3419 // double ytoedge = geom->DetHalfHeight()-std::abs(trkplaney);
3420 // double ztoedge = geom->DetLength()-std::abs(trkplanez);
3421 // std::cout<<"Det W/2,H/2,L "<<geom->DetHalfWidth()<<" "<<geom->DetHalfHeight()<<" "<<geom->DetLength()<<std::endl;
3422 // double cellL = 2*geom->Plane(plane)->Cell(0)->HalfL();
3423 // double cellW = 2*geom->Plane(plane)->Cell(0)->HalfW();
3424 // double cellD = 2*geom->Plane(plane)->Cell(0)->HalfD();
3425 // geo::View_t view = geom->Plane(plane)->View();
3426 /*
3427  if(view == geo::kX){
3428  std::cout<<"x view plane"<<std::endl;
3429  }else if(view == geo::kY){
3430  std::cout<<"y view plane"<<std::endl;
3431  }
3432  std::cout<<"cell L, cellW, cell D "<<cellL<<" "<<cellW<<" "<<cellD<<std::endl;
3433 */
3434  TVector3 toedge(xtoedge, ytoedge, ztoedge);
3435  return toedge;
3436  }
3437 
3438 
3439  bool RecoJMShower::IsFiducial(rb::Prong shower){
3440  int i = shower.ID();
3441  bool isfiducial = true;
3442  for(unsigned int plane=0;plane<fShowerPClusterCol[i].size();plane++){
3443  double planerad = RecoJMShower::Radius(fShowerPClusterCol[i][plane],shower);
3444  TVector3 disttoedge = RecoJMShower::GetTrkPlaneDistToEdge(fShowerPClusterCol[i][plane].Plane(), shower);
3445  if(disttoedge.x()>2*planerad&&disttoedge.y()>2*planerad&&disttoedge.z()>2*planerad)continue;
3446  isfiducial = false;
3447  }
3448  return isfiducial;
3449  }
3450 
3451  double RecoJMShower::GetTrkDistToEdge(rb::Prong shower){
3452  int i = shower.ID();
3453  double mindist = 999999.;
3454  for(unsigned int plane=0;plane<fShowerPClusterCol[i].size();plane++){
3455  TVector3 disttoedge = RecoJMShower::GetTrkPlaneDistToEdge(fShowerPClusterCol[i][plane].Plane(), shower);
3456  double minxyz=999999.;
3457  for(int ii=0;ii<3;ii++){
3458  if(std::abs(disttoedge[ii])<minxyz)minxyz=std::abs(disttoedge[ii]);
3459  }
3460  if(minxyz<mindist)mindist=minxyz;
3461  }
3462  return mindist;
3463  }
3464 
3465 
3466  double RecoJMShower::GetTrkHitPath(unsigned int plane, unsigned int cell, rb::Prong shower){
3467  int i = shower.ID();
3468  Double_t cellXYZ[3]={0,0,0}; // xyz[0] = x, xyz[1] = y, xyz[2] = z
3469  art::ServiceHandle<geo::Geometry> geom;
3470  geom->Plane(plane)->Cell(cell)->GetCenter(cellXYZ);
3471  TVector3 trkdir=shower.Dir();
3472  TVector3 trkstart=shower.Start();
3473  TVector3 trkstop=RecoJMShower::GetShwStop(shower);
3474  double trklength = std::abs((trkstart - trkstop).Mag());
3475  double cellD=2*geom->Plane(plane)->Cell(cell)->HalfD();
3476  double cellL=2*geom->Plane(plane)->Cell(cell)->HalfL();
3477  if(cellL<cellD||cellL<0.0001){
3478  std::cout<<"Geometry error! "<<std::endl;
3479  return 9999;
3480  }
3481  if(std::abs(trkdir.z())< 1e-6||fShowerPClusterCol[i].size()==1){
3482 // std::cout<<"Small Z Dir! "<<std::endl;
3483 // std::cout<<"N Plane = "<<shower.ExtentPlane()<<std::endl;
3484 // std::cout<<"Path = "<<trklength<<std::endl;
3485  return trklength;
3486  }
3487  double trkhitpath=cellD*(sqrt(trkdir.x()*trkdir.x()+trkdir.y()*trkdir.y()+trkdir.z()*trkdir.z())/std::abs(trkdir.z()));
3488  return trkhitpath;
3489  }
3490 
3491 
3492  unsigned int RecoJMShower::ContStartPlane(rb::Prong shower, unsigned int startp, int contp){
3493  int i = shower.ID();
3494  if(int(fShowerPClusterCol[i].size()-startp)<contp)return startp;
3495  unsigned int newstartp = startp;
3496  for(unsigned int plane= startp; plane<fShowerPClusterCol[i].size();plane++){
3497  int starttag = 1;
3498  for(int ic=0;ic<contp;ic++){
3499  if(plane+ic>fShowerPClusterCol[i].size()-1){
3500  starttag=0;
3501  break;
3502  }
3503 // double path = RecoJMShower::GetTrkHitPath(fShowerPClusterCol[i][plane+ic].Plane(), 0, shower);
3504 // double ep = RecoJMShower::Energy(fShowerPClusterCol[i][plane+ic],shower.ID());
3505 // double dedx = ep/path;
3506  double dedx = RecoJMShower::GetPlaneDedx(shower, plane+ic);
3507  if(dedx<0.0005){
3508  starttag=0;
3509  break;
3510  }
3511  }
3512  if(starttag==1){
3513  newstartp = plane;
3514  break;
3515  }
3516  }
3517  return newstartp;
3518  }
3519 
3520 
3521 
3522 //******************************************************
3523 // Energy
3524 //******************************************************
3525 
3526  double RecoJMShower::DepositEnergy(rb::Cluster cluster, int i){
3527  double gev = 0;
3528  for (unsigned int nh=0;nh<cluster.NCell();nh++) {
3529  art::Ptr<rb::CellHit> cellhit = cluster.Cell(nh);
3530  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
3531  double cellEnergy = fHitEWeightMap[key][i];
3532  gev += cellEnergy;
3533  }
3534  return gev;
3535  }
3536 
3537  double RecoJMShower::DepositEnergy(rb::Prong shower){
3538  double gev = 0;
3539  int i = shower.ID();
3540  for (unsigned int nh=0;nh<shower.NCell();nh++) {
3541  art::Ptr<rb::CellHit> cellhit = shower.Cell(nh);
3542  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
3543  double cellEnergy = fHitEWeightMap[key][i];
3544  gev += cellEnergy;
3545  }
3546  return gev;
3547  }
3548 
3549  double RecoJMShower::Energy(rb::Cluster cluster, int i){
3550  double gev = 0;
3551  for (unsigned int nh=0;nh<cluster.NCell();nh++) {
3552  art::Ptr<rb::CellHit> cellhit = cluster.Cell(nh);
3553  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
3554  double cellEnergy = fHitEWeightMap[key][i];
3555  gev += cellEnergy;
3556  }
3557  double ecor = RecoJMShower::ECorrMC(gev);
3558  gev = gev/ecor;
3559  return gev;
3560  }
3561 
3562  double RecoJMShower::Energy(rb::Prong shower){
3563  double gev = 0;
3564  int i = shower.ID();
3565  for (unsigned int nh=0;nh<shower.NCell();nh++) {
3566  art::Ptr<rb::CellHit> cellhit = shower.Cell(nh);
3567  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
3568  double cellEnergy = fHitEWeightMap[key][i];
3569  gev += cellEnergy;
3570  }
3571  double ecor = RecoJMShower::ECorrMC(gev);
3572  double avgx = (shower.Start()[0]+RecoJMShower::GetShwStop(shower)[0])/2.;
3573  double avgy = (shower.Start()[1]+RecoJMShower::GetShwStop(shower)[1])/2.;
3574 // double poscor = (9.9899e-01+7.46983e-05*avgx+7.34449e-05*avgy);//S12-06-17
3575 // double poscor=1.+8.18868e-05*avgx+8.20582e-05*avgy;//without pigtail
3576  double poscor=1.+6.28813e-05*avgx+6.15165e-05*avgy; //with pigtail
3577  if(pCorAbsCell==false) poscor = 1.0+1.04866e-04*avgx+6.40668e-05*avgy;
3578  gev = gev/(ecor*poscor);
3579  //std::cout<<"JMShower Energy calc, depE: "<<showerDepE<<" ecor "<<ecor<<" avgx "<<avgx<<" avgy "<<avgy<<" poscor "<<poscor<<" showerE "<<gev<<std::endl;
3580 // gev = gev/ecor;
3581  return gev;
3582  }
3583 
3584  double RecoJMShower::ECorrMC(double gev){//EM shower energy correction
3585 // double ecor = 5.63616e-01; //without pigtail
3586  double ecor = 5.68142e-01; //with pigtail
3587  if(pCorAbsCell==false)ecor = 5.62343e-01; //no mc cor
3588 // double ecor = 9.05891e-01*(5.53444e-01+9.92889e-02*gev-4.30034e-02*gev*gev-5.55264e-01*exp(-1.61264e+00*gev-1.07062e+01*gev*gev));
3589  return ecor;
3590  }
3591 
3592  double RecoJMShower::CellADCToGeV(double d, geo::View_t view, int detid){
3593  double Att = 4.08050e+03;
3594  if(detid==novadaq::cnv::kFARDET){
3595  if(view==geo::kX){
3596  Att = 4.08574e+03+3.49097e+04*TMath::Exp(-2.67297e-03*d)
3597  -1.62442e+04*TMath::Exp(-3.14363e-02*d-(-3.93281e-05)*d*d-1.49793e-08*d*d*d);
3598  }else if(view == geo::kY){
3599  Att = 4.09347e+03+3.50384e+04*TMath::Exp(-2.68025e-03*d)
3600  -1.74280e+04*TMath::Exp(-3.08694e-02*d-(-3.86870e-05)*d*d-1.48737e-08*d*d*d);
3601  }
3602  }else{
3603  if(view==geo::kX){
3604  Att = 4.08574e+03+3.49097e+04*TMath::Exp(-2.67297e-03*d)
3605  -1.62442e+04*TMath::Exp(-3.14363e-02*d-(-3.93281e-05)*d*d-1.49793e-08*d*d*d);
3606  }else if(view == geo::kY){
3607  Att = 4.09347e+03+3.50384e+04*TMath::Exp(-2.68025e-03*d)
3608  -1.74280e+04*TMath::Exp(-3.08694e-02*d-(-3.86870e-05)*d*d-1.48737e-08*d*d*d);
3609  }
3610  }
3611  return Att;
3612  }
3613 
3614 
3615  double RecoJMShower::CellEDataMC(double gev, double w, geo::View_t view, int detid){//
3616  double ecor = 1; //S12-11-16
3617  if(gev>0&&w>0){
3618  if(detid==novadaq::cnv::kFARDET){
3619  if(view == geo::kX){
3620  ecor = 1.60823+0.000519042*w-2.47044e-07*w*w;
3621  }else if(view == geo::kY){
3622  ecor = 1.3538+0.00051901*w-1.86311e-07*w*w;
3623  }
3624  }else{
3625  if(view == geo::kX){
3626  ecor = 1.60823+0.000519042*w-2.47044e-07*w*w;
3627  }else if(view == geo::kY){
3628  ecor = 1.3538+0.00051901*w-1.86311e-07*w*w;
3629  }
3630  }
3631  }
3632  if(ecor<0.0001)ecor = 1.0;
3633  return ecor;
3634  }
3635 
3636  double RecoJMShower::PlaneDedxDataMC(double gev, int detid){//
3637  double ecor = 1;
3638  if(gev>0){
3639  if(detid==novadaq::cnv::kFARDET){
3640  ecor = 1.0/1.23933;
3641  }else{
3642  ecor = 1.0/1.23933;
3643  }
3644  }
3645  if(ecor<0.0001)ecor = 1.0;
3646  return ecor;
3647  }
3648 
3649  double RecoJMShower::CellTransDataMC(double gev, int detid){//
3650  double ecor = 1;
3651  if(gev>0){
3652  if(detid==novadaq::cnv::kFARDET){
3653  ecor = 1.0/1.336;
3654  }else{
3655  ecor = 1.0/1.336;
3656  }
3657  }
3658  if(ecor<0.0001)ecor = 1.0;
3659  return ecor;
3660  }
3661 
3662  double RecoJMShower::Radius(rb::Prong shower){
3663  int i = shower.ID();
3664  double rad = 0;
3665  double gev = 0;
3666  for (unsigned int nh=0;nh<shower.NCell();nh++) {
3667  art::Ptr<rb::CellHit> cellhit = shower.Cell(nh);
3668  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
3669  double cellEnergy = fHitEWeightMap[key][i];
3670 // double celldist = (RecoJMShower::GetTrkHitPos(cellhit->Plane(),cellhit->Cell(),shower.Dir(),shower.Start(),RecoJMShower::GetShwStop(shower))).z();
3671  double celldist= RecoJMShower::GetCellDistToTrk(cellhit->Plane(),cellhit->Cell(),shower.Start(),RecoJMShower::GetShwStop(shower));
3672  rad += fHitEWeightMap[key][i]*celldist;
3673  gev += cellEnergy;
3674  }
3675  if(gev>0.001)rad = rad/gev;
3676  return rad;
3677  }
3678 
3679  double RecoJMShower::Radius(rb::Cluster cluster, rb::Prong shower){
3680  int i = shower.ID();
3681  double rad = 0;
3682  double gev = 0;
3683  for (unsigned int nh=0;nh<cluster.NCell();nh++) {
3684  art::Ptr<rb::CellHit> cellhit = cluster.Cell(nh);
3685  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
3686  double cellEnergy = fHitEWeightMap[key][i];
3687 // double celldist = (RecoJMShower::GetTrkHitPos(cellhit->Plane(),cellhit->Cell(),shower.Dir(),shower.Start(),RecoJMShower::GetShwStop(shower))).z();
3688  double celldist= RecoJMShower::GetCellDistToTrk(cellhit->Plane(),cellhit->Cell(),shower.Start(),RecoJMShower::GetShwStop(shower));
3689  rad += fHitEWeightMap[key][i]*celldist;
3690  gev += cellEnergy;
3691  }
3692  if(gev>0.001)rad = rad/gev;
3693  return rad;
3694  }
3695 
3696  double RecoJMShower::RadiusV(rb::Cluster cluster){
3697  double rad = 0;
3698  double gev = 0;
3699  TVector3 clustercentroid = RecoJMShower::GetCentroid(cluster);
3700  for (unsigned int nh=0;nh<cluster.NCell();nh++) {
3701  art::Ptr<rb::CellHit> cellhit = cluster.Cell(nh);
3702  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
3703  double cellEnergy = cellhit->PE();
3704  TVector3 celldist = RecoJMShower::GetCellDistToPoint(cellhit->Plane(), cellhit->Cell(),clustercentroid);
3705  rad += cellEnergy*sqrt(celldist[0]*celldist[0]+celldist[1]*celldist[1]+celldist[2]*celldist[2]);
3706  gev += cellEnergy;
3707  }
3708  if(gev>0.001)rad = rad/gev;
3709  return rad;
3710  }
3711 
3712  TVector3 RecoJMShower::GetCentroid(rb::Cluster cluster){
3713  double xadc = 0;
3714  double yadc = 0;
3715  double zadc = 0;
3716  double xadccellx = 0;
3717  double yadccelly = 0;
3718  double zadccellz = 0;
3719  double cx = 0;
3720  double cy = 0;
3721  double cz = 0;
3722 
3723  for(unsigned int nh=0; nh<cluster.NXCell();nh++){
3724  art::Ptr<rb::CellHit> cellhit = cluster.XCell(nh);
3725  Double_t cellXYZ[3]={0,0,0};
3726  art::ServiceHandle<geo::Geometry> geom;
3727  geom->Plane(cellhit->Plane())->Cell(cellhit->Cell())->GetCenter(cellXYZ);
3728  double cellADC = double(cellhit->ADC());
3729  xadc += cellADC;
3730  zadc += cellADC;
3731  xadccellx += cellADC*cellXYZ[0];
3732  zadccellz += cellADC*cellXYZ[2];
3733  }
3734 
3735  for(unsigned int nh=0; nh<cluster.NYCell();nh++){
3736  art::Ptr<rb::CellHit> cellhit = cluster.YCell(nh);
3737  Double_t cellXYZ[3]={0,0,0};
3738  art::ServiceHandle<geo::Geometry> geom;
3739  geom->Plane(cellhit->Plane())->Cell(cellhit->Cell())->GetCenter(cellXYZ);
3740  double cellADC = double(cellhit->ADC());
3741  yadc += cellADC;
3742  zadc += cellADC;
3743  yadccelly += cellADC*cellXYZ[1];
3744  zadccellz += cellADC*cellXYZ[2];
3745  }
3746  if(xadc>0)cx = xadccellx/xadc;
3747  if(yadc>0)cy = yadccelly/yadc;
3748  if(zadc>0)cz = zadccellz/zadc;
3749  TVector3 centroid(cx, cy, cz);
3750  return centroid;
3751  }
3752 
3753 
3754  double RecoJMShower::Gap(rb::Prong shower, TVector3 vertex){
3755  TVector3 showerStart=shower.Start();
3756  double gap = (vertex-showerStart).Mag();
3757  return gap;
3758  }
3759 
3760  double RecoJMShower::GetPlaneDedx(rb::Prong shower, unsigned int plane){
3761  int i = shower.ID();
3762  if(fShowerPClusterCol[i].size()<1)
3763  {
3764  std::cout<<"JMShower: Nothing in plane "<<i<<std::endl;
3765  return 9999;
3766  }
3767  double path = RecoJMShower::GetTrkHitPath(fShowerPClusterCol[i][plane].Plane(), 0, shower);
3768  double ep = RecoJMShower::Energy(fShowerPClusterCol[i][plane], shower.ID());
3769 // TVector3 hitpos = RecoJMShower::GetTrkPlanePos(fShowerPClusterCol[i][plane].Plane(), shower);
3770 // double poscor = (0.927377+6.46394e-05*hitpos[0])*(0.993889+3.83504e-05*hitpos[1]);/from plane
3771 // double poscor = (9.9899+7.46983e-05*avgx+7.34449e-05*avgy);//from shower
3772 // double dedx = (ep/poscor)/path;
3773  double dedx = ep/path;
3774  if(isrealdata==true&&pCorDataDedx)dedx=dedx/RecoJMShower::PlaneDedxDataMC(dedx, fDetid);
3775  //std::cout<<"shower: "<<i<<" plane idx: "<<plane<<" plane: "<<fShowerPClusterCol[i][plane].Plane()<<" path: "<<path<<" ep: "<<ep<<" dedx: "<<dedx<<std::endl;
3776  return dedx;
3777  }
3778 
3779 
3780  double RecoJMShower::GetPlaneDedxProb(rb::Prong shower, int ireg, int igev, unsigned int plane, double dedx, int type){
3781 // int i = shower.ID();
3782  double p = 0;
3783 
3784  if(int(jmshower::kElectronSG)==type){
3785  if(ireg<0||ireg>3)return 1;
3786  if(igev>=41||igev<0)return 1;
3787  if(plane>=200)return 1;
3788  }else{
3789  if(ireg<0||ireg>3)return 1;
3790  if(igev>=11||igev<0)return 1;
3791  if(plane>=200)return 1;
3792  }
3793 
3794  if(int(jmshower::kElectron)==type){
3795  double ntot = double(fHtPlaneDedxE[ireg][igev][plane]->GetEntries());
3796  double n = double(fHtPlaneDedxE[ireg][igev][plane]->Interpolate(dedx));
3797  double nbin = fHtPlaneDedxE[ireg][igev][plane]->GetNbinsX();
3798  if(n<0.0001)n=0.0001;
3799  if(ntot<1){
3800  p = 0;
3801  }else{
3802  p = n*nbin/ntot;
3803  }
3804  //std::cout<<"type "<<type<<" "<<ireg<<" "<<igev<<" "<<plane<<" "<<dedx<<" "<<ntot<<" "<<n<<" "<<nbin<<" "<<p<<std::endl;
3805  }
3806  if(int(jmshower::kGamma)==type){
3807  double ntot = double(fHtPlaneDedxG[ireg][igev][plane]->GetEntries());
3808  double n = double(fHtPlaneDedxG[ireg][igev][plane]->Interpolate(dedx));
3809  double nbin = fHtPlaneDedxG[ireg][igev][plane]->GetNbinsX();
3810  if(n<0.0001)n=0.0001;
3811  if(ntot<1){
3812  p = 0;
3813  }else{
3814  p = n*nbin/ntot;
3815  }
3816  }
3817  if(int(jmshower::kMuon)==type){
3818  double ntot = double(fHtPlaneDedxMu[ireg][igev][plane]->GetEntries());
3819  double n = double(fHtPlaneDedxMu[ireg][igev][plane]->Interpolate(dedx));
3820  double nbin = fHtPlaneDedxMu[ireg][igev][plane]->GetNbinsX();
3821  if(n<0.0001)n=0.0001;
3822  if(ntot<1){
3823  p = 0;
3824  }else{
3825  p = n*nbin/ntot;
3826  }
3827  }
3828  if(int(jmshower::kPi0)==type){
3829  double ntot = double(fHtPlaneDedxPi0[ireg][igev][plane]->GetEntries());
3830  double n = double(fHtPlaneDedxPi0[ireg][igev][plane]->Interpolate(dedx));
3831  double nbin = fHtPlaneDedxPi0[ireg][igev][plane]->GetNbinsX();
3832  if(n<0.0001)n=0.0001;
3833  if(ntot<1){
3834  p = 0;
3835  }else{
3836  p = n*nbin/ntot;
3837  }
3838  }
3839 /*
3840  if(int(jmshower::kHad)==type){
3841  double ntot = double(fHtPlaneDedxHad[ireg][igev][plane]->GetEntries());
3842  double n = double(fHtPlaneDedxHad[ireg][igev][plane]->Interpolate(dedx));
3843  double nbin = fHtPlaneDedxHad[ireg][igev][plane]->GetNbinsX();
3844  if(n<0.0001)n=0.0001;
3845  if(ntot<1){
3846  p = 0;
3847  }else{
3848  p = n*nbin/ntot;
3849  }
3850  }
3851 */
3852 
3853  if(int(jmshower::kProton)==type){
3854  double ntot = double(fHtPlaneDedxP[ireg][igev][plane]->GetEntries());
3855  double n = double(fHtPlaneDedxP[ireg][igev][plane]->Interpolate(dedx));
3856  double nbin = fHtPlaneDedxP[ireg][igev][plane]->GetNbinsX();
3857  if(n<0.0001)n=0.0001;
3858  if(ntot<1){
3859  p = 0;
3860  }else{
3861  p = n*nbin/ntot;
3862  }
3863  }
3864 
3865  if(int(jmshower::kNeutron)==type){
3866  double ntot = double(fHtPlaneDedxN[ireg][igev][plane]->GetEntries());
3867  double n = double(fHtPlaneDedxN[ireg][igev][plane]->Interpolate(dedx));
3868  double nbin = fHtPlaneDedxN[ireg][igev][plane]->GetNbinsX();
3869  if(n<0.0001)n=0.0001;
3870  if(ntot<1){
3871  p = 0;
3872  }else{
3873  p = n*nbin/ntot;
3874  }
3875  }
3876 
3877  if(int(jmshower::kPion)==type){
3878  double ntot = double(fHtPlaneDedxPi[ireg][igev][plane]->GetEntries());
3879  double n = double(fHtPlaneDedxPi[ireg][igev][plane]->Interpolate(dedx));
3880  double nbin = fHtPlaneDedxPi[ireg][igev][plane]->GetNbinsX();
3881  if(n<0.0001)n=0.0001;
3882  if(ntot<1){
3883  p = 0;
3884  }else{
3885  p = n*nbin/ntot;
3886  }
3887  }
3888  if(int(jmshower::kElectronQE)==type){
3889  double ntot = double(fHtPlaneDedxEqe[ireg][igev][plane]->GetEntries());
3890  double n = double(fHtPlaneDedxEqe[ireg][igev][plane]->Interpolate(dedx));
3891  double nbin = fHtPlaneDedxEqe[ireg][igev][plane]->GetNbinsX();
3892  if(n<0.0001)n=0.0001;
3893  if(ntot<1){
3894  p = 0;
3895  }else{
3896  p = n*nbin/ntot;
3897  }
3898  }
3899 
3900  if(int(jmshower::kElectronRES)==type){
3901  double ntot = double(fHtPlaneDedxEres[ireg][igev][plane]->GetEntries());
3902  double n = double(fHtPlaneDedxEres[ireg][igev][plane]->Interpolate(dedx));
3903  double nbin = fHtPlaneDedxEres[ireg][igev][plane]->GetNbinsX();
3904  if(n<0.0001)n=0.0001;
3905  if(ntot<1){
3906  p = 0;
3907  }else{
3908  p = n*nbin/ntot;
3909  }
3910  }
3911 
3912  if(int(jmshower::kElectronDIS)==type){
3913  double ntot = double(fHtPlaneDedxEdis[ireg][igev][plane]->GetEntries());
3914  double n = double(fHtPlaneDedxEdis[ireg][igev][plane]->Interpolate(dedx));
3915  double nbin = fHtPlaneDedxEdis[ireg][igev][plane]->GetNbinsX();
3916  if(n<0.0001)n=0.0001;
3917  if(ntot<1){
3918  p = 0;
3919  }else{
3920  p = n*nbin/ntot;
3921  }
3922  }
3923 
3924  if(int(jmshower::kElectronCOH)==type){
3925  double ntot = double(fHtPlaneDedxEcoh[ireg][igev][plane]->GetEntries());
3926  double n = double(fHtPlaneDedxEcoh[ireg][igev][plane]->Interpolate(dedx));
3927  double nbin = fHtPlaneDedxEcoh[ireg][igev][plane]->GetNbinsX();
3928  if(n<0.0001)n=0.0001;
3929  if(ntot<1){
3930  p = 0;
3931  }else{
3932  p = n*nbin/ntot;
3933  }
3934  }
3935 
3936  if(int(jmshower::kElectronSG)==type){
3937  double ntot = double(fHtPlaneDedxEsg[ireg][igev][plane]->GetEntries());
3938  double n = double(fHtPlaneDedxEsg[ireg][igev][plane]->Interpolate(dedx));
3939  double nbin = fHtPlaneDedxEsg[ireg][igev][plane]->GetNbinsX();
3940  if(n<0.0001)n=0.0001;
3941  if(ntot<1){
3942  p = 0;
3943  }else{
3944  p = n*nbin/ntot;
3945  }
3946  }
3947  return p;
3948  }
3949 
3950 /*
3951  double RecoJMShower::GetPlaneDedxProb(rb::Prong shower, int igev, unsigned int plane, double dedx, unsigned int nmip, const int type){
3952 // int i = shower.ID();
3953  double p = 0;
3954  if(igev>=11||igev<0)return 1;
3955  if(plane>=200)return 1;
3956 // double dedx1 = RecoJMShower::GetPlaneDedx(shower, plane+1);
3957  if(int(jmshower::kElectron)==type){
3958  if(plane<nmip){
3959  double ntot = double(fHtmipPlaneDedxE[igev][plane]->GetEntries());
3960  double ntotcomb = double(fHtmipPlaneDedxE[igev][plane]->GetEntries()+fHtshowerPlaneDedxE[igev][plane]->GetEntries());
3961  double n = double(fHtmipPlaneDedxE[igev][plane]->Interpolate(dedx));
3962  double nbin = fHtmipPlaneDedxE[igev][plane]->GetNbinsX();
3963 // double nbincomb = fHtmipPlaneDedxE[igev][plane]->GetNbinsX()+fHtshowerPlaneDedxE[igev][plane]->GetNbinsX();
3964  if(n<0.0001)n=0.0001;
3965  if(ntotcomb<100){
3966  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
3967  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
3968  }else{
3969  p = (n*nbin/ntot)*(ntot/ntotcomb);
3970  }
3971  }
3972  if(plane>=nmip){
3973 // double ntot = double(fHtshowerPlaneDedxE[igev][plane]->GetEntries());
3974  double ntotcomb = double(fHtshowerPlaneDedxE[igev][plane]->GetEntries()+fHtmipPlaneDedxE[igev][plane]->GetEntries());
3975  double n = double(fHtshowerPlaneDedxE[igev][plane]->Interpolate(dedx));
3976  double nbin = fHtshowerPlaneDedxE[igev][plane]->GetNbinsX();
3977 // double nbincomb = fHtshowerPlaneDedxE[igev][plane]->GetNbinsX()+fHtmipPlaneDedxE[igev][plane]->GetNbinsX();
3978  if(n<0.0001)n=0.0001;
3979  if(ntotcomb<100){
3980  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
3981  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
3982  }else{
3983  p = n*nbin/ntotcomb;
3984  }
3985  }
3986  }
3987 
3988 
3989  if(int(jmshower::kGamma)==type){
3990  if(plane<nmip){
3991 // double ntot = double(fHtmipPlaneDedxG[igev][plane]->GetEntries());
3992  double ntotcomb = double(fHtmipPlaneDedxG[igev][plane]->GetEntries()+fHtshowerPlaneDedxG[igev][plane]->GetEntries());
3993  double n = double(fHtmipPlaneDedxG[igev][plane]->Interpolate(dedx));
3994  double nbin = double(fHtmipPlaneDedxG[igev][plane]->GetNbinsX());
3995 // double nbincomb = double(fHtmipPlaneDedxG[igev][plane]->GetNbinsX()+fHtshowerPlaneDedxG[igev][plane]->GetNbinsX());
3996  if(n<0.0001)n=0.0001;
3997  if(ntotcomb<100){
3998  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
3999  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4000  }else{
4001  p = n*nbin/ntotcomb;
4002  }
4003  }
4004  if(plane>=nmip){
4005 // double ntot = double(fHtshowerPlaneDedxG[igev][plane]->GetEntries());
4006  double ntotcomb = double(fHtshowerPlaneDedxG[igev][plane]->GetEntries()+fHtmipPlaneDedxG[igev][plane]->GetEntries());
4007  double n = double(fHtshowerPlaneDedxG[igev][plane]->Interpolate(dedx));
4008  double nbin = double(fHtshowerPlaneDedxG[igev][plane]->GetNbinsX());
4009 // double nbincomb = double(fHtshowerPlaneDedxG[igev][plane]->GetNbinsX()+fHtmipPlaneDedxG[igev][plane]->GetNbinsX());
4010  if(n<0.0001)n=0.0001;
4011  if(ntotcomb<100){
4012  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4013  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4014  }else{
4015  p = n*nbin/ntotcomb;
4016  }
4017  }
4018  }
4019 
4020 
4021  if(int(jmshower::kMuon)==type){
4022  double ntot = double(fHtPlaneDedxMu[igev][plane]->GetEntries());
4023  double n = double(fHtPlaneDedxMu[igev][plane]->Interpolate(dedx));
4024  double nbin = fHtPlaneDedxMu[igev][plane]->GetNbinsX();
4025  if(n<0.0001)n=0.0001;
4026  if(ntot<1&&igev>0&&igev<10){
4027  p = 0;
4028  }else{
4029  p = n*nbin/ntot;
4030  }
4031  }
4032 
4033  if(int(jmshower::kPi0)==type){
4034  if(plane<nmip){
4035 // double ntot = double(fHtmipPlaneDedxPi0[igev][plane]->GetEntries());
4036  double ntotcomb = double(fHtmipPlaneDedxPi0[igev][plane]->GetEntries()+fHtshowerPlaneDedxPi0[igev][plane]->GetEntries());
4037  double n = double(fHtmipPlaneDedxPi0[igev][plane]->Interpolate(dedx));
4038  double nbin = double(fHtmipPlaneDedxPi0[igev][plane]->GetNbinsX());
4039 // double nbincomb = double(fHtmipPlaneDedxPi0[igev][plane]->GetNbinsX()+fHtshowerPlaneDedxPi0[igev][plane]->GetNbinsX());
4040  if(n<0.0001)n=0.0001;
4041  if(ntotcomb<100){
4042  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4043  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4044  }else{
4045  p = n*nbin/ntotcomb;
4046  }
4047  }
4048  if(plane>=nmip){
4049 // double ntot = double(fHtshowerPlaneDedxPi0[igev][plane]->GetEntries());
4050  double ntotcomb = double(fHtshowerPlaneDedxPi0[igev][plane]->GetEntries()+fHtmipPlaneDedxPi0[igev][plane]->GetEntries());
4051  double n = double(fHtshowerPlaneDedxPi0[igev][plane]->Interpolate(dedx));
4052  double nbin = double(fHtshowerPlaneDedxPi0[igev][plane]->GetNbinsX());
4053 // double nbincomb = double(fHtshowerPlaneDedxPi0[igev][plane]->GetNbinsX()+fHtmipPlaneDedxPi0[igev][plane]->GetNbinsX());
4054  if(n<0.0001)n=0.0001;
4055  if(ntotcomb<100){
4056  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4057  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4058  }else{
4059  p = n*nbin/ntotcomb;
4060  }
4061  }
4062  }
4063 
4064 
4065  if(int(jmshower::kHad)==type){
4066  if(plane<nmip){
4067 // double ntot = double(fHtmipPlaneDedxHad[igev][plane]->GetEntries());
4068  double ntotcomb = double(fHtmipPlaneDedxHad[igev][plane]->GetEntries()+fHtshowerPlaneDedxHad[igev][plane]->GetEntries());
4069  double n = double(fHtmipPlaneDedxHad[igev][plane]->Interpolate(dedx));
4070  double nbin = double(fHtmipPlaneDedxHad[igev][plane]->GetNbinsX());
4071 // double nbincomb = double(fHtmipPlaneDedxHad[igev][plane]->GetNbinsX()+fHtshowerPlaneDedxHad[igev][plane]->GetNbinsX());
4072  if(n<0.0001)n=0.0001;
4073  if(ntotcomb<100){
4074  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4075  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4076  }else{
4077  p = n*nbin/ntotcomb;
4078  }
4079  }
4080  if(plane>=nmip){
4081 // double ntot = double(fHtshowerPlaneDedxHad[igev][plane]->GetEntries());
4082  double ntotcomb = double(fHtshowerPlaneDedxHad[igev][plane]->GetEntries()+fHtmipPlaneDedxHad[igev][plane]->GetEntries());
4083  double n = double(fHtshowerPlaneDedxHad[igev][plane]->Interpolate(dedx));
4084  double nbin = double(fHtshowerPlaneDedxHad[igev][plane]->GetNbinsX());
4085 // double nbincomb = double(fHtshowerPlaneDedxHad[igev][plane]->GetNbinsX()+fHtmipPlaneDedxHad[igev][plane]->GetNbinsX());
4086  if(n<0.0001)n=0.0001;
4087  if(ntotcomb<100){
4088  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4089  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4090  }else{
4091  p = n*nbin/ntotcomb;
4092  }
4093  }
4094  }
4095 
4096 
4097  if(int(jmshower::kProton)==type){
4098  if(plane<nmip){
4099 // double ntot = double(fHtmipPlaneDedxP[igev][plane]->GetEntries());
4100  double ntotcomb = double(fHtmipPlaneDedxP[igev][plane]->GetEntries()+fHtshowerPlaneDedxP[igev][plane]->GetEntries());
4101  double n = double(fHtmipPlaneDedxP[igev][plane]->Interpolate(dedx));
4102  double nbin = double(fHtmipPlaneDedxP[igev][plane]->GetNbinsX());
4103 // double nbincomb = double(fHtmipPlaneDedxP[igev][plane]->GetNbinsX()+fHtshowerPlaneDedxP[igev][plane]->GetNbinsX());
4104  if(n<0.0001)n=0.0001;
4105  if(ntotcomb<100){
4106  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4107  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4108  }else{
4109  p = n*nbin/ntotcomb;
4110  }
4111  }
4112  if(plane>=nmip){
4113 // double ntot = double(fHtshowerPlaneDedxP[igev][plane]->GetEntries());
4114  double ntotcomb = double(fHtshowerPlaneDedxP[igev][plane]->GetEntries()+fHtmipPlaneDedxP[igev][plane]->GetEntries());
4115  double n = double(fHtshowerPlaneDedxP[igev][plane]->Interpolate(dedx));
4116  double nbin = double(fHtshowerPlaneDedxP[igev][plane]->GetNbinsX());
4117 // double nbincomb = double(fHtshowerPlaneDedxP[igev][plane]->GetNbinsX()+fHtmipPlaneDedxP[igev][plane]->GetNbinsX());
4118  if(n<0.0001)n=0.0001;
4119  if(ntotcomb<100){
4120  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4121  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4122  }else{
4123  p = n*nbin/ntotcomb;
4124  }
4125  }
4126  }
4127 
4128 
4129 
4130  if(int(jmshower::kNeutron)==type){
4131  if(plane<nmip){
4132 // double ntot = double(fHtmipPlaneDedxN[igev][plane]->GetEntries());
4133  double ntotcomb = double(fHtmipPlaneDedxN[igev][plane]->GetEntries()+fHtshowerPlaneDedxN[igev][plane]->GetEntries());
4134  double n = double(fHtmipPlaneDedxN[igev][plane]->Interpolate(dedx));
4135  double nbin = double(fHtmipPlaneDedxN[igev][plane]->GetNbinsX());
4136 // double nbincomb = double(fHtmipPlaneDedxN[igev][plane]->GetNbinsX()+fHtshowerPlaneDedxN[igev][plane]->GetNbinsX());
4137  if(n<0.0001)n=0.0001;
4138  if(ntotcomb<100){
4139  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4140  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4141  }else{
4142  p = n*nbin/ntotcomb;
4143  }
4144  }
4145  if(plane>=nmip){
4146 // double ntot = double(fHtshowerPlaneDedxN[igev][plane]->GetEntries());
4147  double ntotcomb = double(fHtshowerPlaneDedxN[igev][plane]->GetEntries()+fHtmipPlaneDedxN[igev][plane]->GetEntries());
4148  double n = double(fHtshowerPlaneDedxN[igev][plane]->Interpolate(dedx));
4149  double nbin = double(fHtshowerPlaneDedxN[igev][plane]->GetNbinsX());
4150 // double nbincomb = double(fHtshowerPlaneDedxN[igev][plane]->GetNbinsX()+fHtmipPlaneDedxN[igev][plane]->GetNbinsX());
4151  if(n<0.0001)n=0.0001;
4152  if(ntotcomb<100){
4153  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4154  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4155  }else{
4156  p = n*nbin/ntotcomb;
4157  }
4158  }
4159  }
4160 
4161  if(int(jmshower::kPion)==type){
4162  if(plane<nmip){
4163 // double ntot = double(fHtmipPlaneDedxPi[igev][plane]->GetEntries());
4164  double ntotcomb = double(fHtmipPlaneDedxPi[igev][plane]->GetEntries()+fHtshowerPlaneDedxPi[igev][plane]->GetEntries());
4165  double n = double(fHtmipPlaneDedxPi[igev][plane]->Interpolate(dedx));
4166  double nbin = double(fHtmipPlaneDedxPi[igev][plane]->GetNbinsX());
4167 // double nbincomb = double(fHtmipPlaneDedxPi[igev][plane]->GetNbinsX()+fHtshowerPlaneDedxPi[igev][plane]->GetNbinsX());
4168  if(n<0.0001)n=0.0001;
4169  if(ntotcomb<100){
4170  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4171  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4172  }else{
4173  p = n*nbin/ntotcomb;
4174  }
4175  }
4176  if(plane>=nmip){
4177 // double ntot = double(fHtshowerPlaneDedxPi[igev][plane]->GetEntries());
4178  double ntotcomb = double(fHtshowerPlaneDedxPi[igev][plane]->GetEntries()+fHtmipPlaneDedxPi[igev][plane]->GetEntries());
4179  double n = double(fHtshowerPlaneDedxPi[igev][plane]->Interpolate(dedx));
4180  double nbin = double(fHtshowerPlaneDedxPi[igev][plane]->GetNbinsX());
4181 // double nbincomb = double(fHtshowerPlaneDedxPi[igev][plane]->GetNbinsX()+fHtmipPlaneDedxPi[igev][plane]->GetNbinsX());
4182  if(n<0.0001)n=0.0001;
4183  if(ntotcomb<100){
4184  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4185  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4186  }else{
4187  p = n*nbin/ntotcomb;
4188  }
4189  }
4190  }
4191 
4192 
4193 
4194  if(int(jmshower::kElectronQE)==type){
4195  if(plane<nmip){
4196  double ntotcomb = double(fHtmipPlaneDedxEqe[igev][plane]->GetEntries()+fHtshowerPlaneDedxEqe[igev][plane]->GetEntries());
4197  double n = double(fHtmipPlaneDedxEqe[igev][plane]->Interpolate(dedx));
4198  double nbin = double(fHtmipPlaneDedxEqe[igev][plane]->GetNbinsX());
4199  if(n<0.0001)n=0.0001;
4200  if(ntotcomb<100){
4201  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4202  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4203  }else{
4204  p = n*nbin/ntotcomb;
4205  }
4206  }
4207  if(plane>=nmip){
4208  double ntotcomb = double(fHtshowerPlaneDedxEqe[igev][plane]->GetEntries()+fHtmipPlaneDedxEqe[igev][plane]->GetEntries());
4209  double n = double(fHtshowerPlaneDedxEqe[igev][plane]->Interpolate(dedx));
4210  double nbin = double(fHtshowerPlaneDedxEqe[igev][plane]->GetNbinsX());
4211  if(n<0.0001)n=0.0001;
4212  if(ntotcomb<100){
4213  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4214  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4215  }else{
4216  p = n*nbin/ntotcomb;
4217  }
4218  }
4219  }
4220 
4221 
4222  if(int(jmshower::kElectronRES)==type){
4223  if(plane<nmip){
4224  double ntotcomb = double(fHtmipPlaneDedxEres[igev][plane]->GetEntries()+fHtshowerPlaneDedxEres[igev][plane]->GetEntries());
4225  double n = double(fHtmipPlaneDedxEres[igev][plane]->Interpolate(dedx));
4226  double nbin = double(fHtmipPlaneDedxEres[igev][plane]->GetNbinsX());
4227  if(n<0.0001)n=0.0001;
4228  if(ntotcomb<100){
4229  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4230  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4231  }else{
4232  p = n*nbin/ntotcomb;
4233  }
4234  }
4235  if(plane>=nmip){
4236  double ntotcomb = double(fHtshowerPlaneDedxEres[igev][plane]->GetEntries()+fHtmipPlaneDedxEres[igev][plane]->GetEntries());
4237  double n = double(fHtshowerPlaneDedxEres[igev][plane]->Interpolate(dedx));
4238  double nbin = double(fHtshowerPlaneDedxEres[igev][plane]->GetNbinsX());
4239  if(n<0.0001)n=0.0001;
4240  if(ntotcomb<100){
4241  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4242  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4243  }else{
4244  p = n*nbin/ntotcomb;
4245  }
4246  }
4247  }
4248 
4249 
4250  if(int(jmshower::kElectronDIS)==type){
4251  if(plane<nmip){
4252  double ntotcomb = double(fHtmipPlaneDedxEdis[igev][plane]->GetEntries()+fHtshowerPlaneDedxEdis[igev][plane]->GetEntries());
4253  double n = double(fHtmipPlaneDedxEdis[igev][plane]->Interpolate(dedx));
4254  double nbin = double(fHtmipPlaneDedxEdis[igev][plane]->GetNbinsX());
4255  if(n<0.0001)n=0.0001;
4256  if(ntotcomb<100){
4257  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4258  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4259  }else{
4260  p = n*nbin/ntotcomb;
4261  }
4262  }
4263  if(plane>=nmip){
4264  double ntotcomb = double(fHtshowerPlaneDedxEdis[igev][plane]->GetEntries()+fHtmipPlaneDedxEdis[igev][plane]->GetEntries());
4265  double n = double(fHtshowerPlaneDedxEdis[igev][plane]->Interpolate(dedx));
4266  double nbin = double(fHtshowerPlaneDedxEdis[igev][plane]->GetNbinsX());
4267  if(n<0.0001)n=0.0001;
4268  if(ntotcomb<100){
4269  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4270  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4271  }else{
4272  p = n*nbin/ntotcomb;
4273  }
4274  }
4275  }
4276 
4277 
4278  if(int(jmshower::kElectronCOH)==type){
4279  if(plane<nmip){
4280  double ntotcomb = double(fHtmipPlaneDedxEcoh[igev][plane]->GetEntries()+fHtshowerPlaneDedxEcoh[igev][plane]->GetEntries());
4281  double n = double(fHtmipPlaneDedxEcoh[igev][plane]->Interpolate(dedx));
4282  double nbin = double(fHtmipPlaneDedxEcoh[igev][plane]->GetNbinsX());
4283  if(n<0.0001)n=0.0001;
4284  if(ntotcomb<100){
4285  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4286  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4287  }else{
4288  p = n*nbin/ntotcomb;
4289  }
4290  }
4291  if(plane>=nmip){
4292  double ntotcomb = double(fHtshowerPlaneDedxEcoh[igev][plane]->GetEntries()+fHtmipPlaneDedxEcoh[igev][plane]->GetEntries());
4293  double n = double(fHtshowerPlaneDedxEcoh[igev][plane]->Interpolate(dedx));
4294  double nbin = double(fHtshowerPlaneDedxEcoh[igev][plane]->GetNbinsX());
4295  if(n<0.0001)n=0.0001;
4296  if(ntotcomb<100){
4297  if(igev>1)p = RecoJMShower::GetPlaneDedxProb(shower, igev-1, plane, dedx, nmip, type);
4298  if(igev<=1)p = RecoJMShower::GetPlaneDedxProb(shower, igev+1, plane, dedx, nmip, type);
4299  }else{
4300  p = n*nbin/ntotcomb;
4301  }
4302  }
4303  }
4304 
4305 
4306  return p;
4307  }
4308 */
4309 
4310  double RecoJMShower::GetInterPlaneDedxProb(rb::Prong shower, unsigned int plane, double dedx, const int type, TVector3 pos, int iReg, int igev0, int igev1, double e0, double e1){
4311 
4312  double probPXPY = 0;
4313  bool interPos = false;
4314 
4315 // int i = shower.ID();
4316  art::ServiceHandle<geo::Geometry> geom;
4317 // TVector3 pos = RecoJMShower::GetTrkPlanePos(fShowerPClusterCol[i][plane].Plane(), shower);
4318 // std::cout<<"in getinterplane pos X,Y "<<pos[0]<<" "<<pos[1]<<std::endl;
4319 
4320 // TVector3 pos((shower.Start()[0]+RecoJMShower::GetShwStop(shower)[0])/2.0,(shower.Start()[1]+RecoJMShower::GetShwStop(shower)[1])/2.0,(shower.Start()[2]+RecoJMShower::GetShwStop(shower)[2])/2.0);
4321 
4322 
4323  double cos = shower.Dir()[2];
4324 
4325  if(interPos==true){
4326  double X0 = -geom->DetHalfWidth()/2.0;
4327  double X1 = geom->DetHalfWidth()/2.0;
4328  double Y0 = -geom->DetHalfHeight()/2.0;
4329  double Y1 = geom->DetHalfHeight()/2.0;
4330 
4331  double probreg[4]={1,1,1,1};
4332  for(int ireg=0;ireg<4;ireg++){
4333  unsigned int plane0 = int(int(plane)/std::abs(cos));
4334  unsigned int plane1 = int(int(plane)/std::abs(cos)) + 1;
4335  double r0 = double(plane/std::abs(cos))-double(plane0);
4336  double r1 = 1-r0;
4337  double probe0r0 = 1;
4338  double probe0r1 = 1;
4339  double probe1r0 = 1;
4340  double probe1r1 = 1;
4341 
4342  probe0r0 = RecoJMShower::GetPlaneDedxProb(shower, ireg, igev0, plane0, dedx, type);
4343  probe0r1 = RecoJMShower::GetPlaneDedxProb(shower, ireg, igev0, plane1, dedx, type);
4344  probe1r0 = RecoJMShower::GetPlaneDedxProb(shower, ireg, igev1, plane0, dedx, type);
4345  probe1r1 = RecoJMShower::GetPlaneDedxProb(shower, ireg, igev1, plane1, dedx, type);
4346 // std::cout<<"ireg, prob probe0r0 "<<ireg<<" "<<probe0r0<<std::endl;
4347 // std::cout<<"ireg, prob probe0r1 "<<ireg<<" "<<probe0r1<<std::endl;
4348 // std::cout<<"ireg, prob probe1r0 "<<ireg<<" "<<probe1r0<<std::endl;
4349 // std::cout<<"ireg, prob probe1r0 "<<ireg<<" "<<probe1r0<<std::endl;
4350  double probe0 = (r1*probe0r0+r0*probe0r1);
4351  double probe1 = (r1*probe1r0+r0*probe1r1);
4352  if(igev0!=igev1)probreg[ireg] = (e1*probe0+e0*probe1)/(e0+e1);
4353  if(igev0==igev1)probreg[ireg] = probe0;
4354 // std::cout<<"ireg, probreg "<<ireg<<" "<<probreg[ireg]<<std::endl;
4355  }
4356  double PX = pos[0];
4357  double PY = pos[1];
4358  double probX1Y1 = probreg[0];
4359  double probX1Y0 = probreg[1];
4360  double probX0Y1 = probreg[2];
4361  double probX0Y0 = probreg[3];
4362 // std::cout<<"prob X1Y1, X1Y0, X0Y1, X0Y0 "<<probX1Y1<<" "<<probX1Y0<<" "<<probX0Y1<<" "<<probX0Y0<<std::endl;
4363  double probX1PY = ((PY-Y0)/(Y1-Y0))*(probX1Y1-probX1Y0)+probX1Y0;
4364  double probX0PY = ((PY-Y0)/(Y1-Y0))*(probX0Y1-probX0Y0)+probX0Y0;
4365  probPXPY = ((PX-X0)/(X1-X0))*(probX1PY-probX0PY)+probX0PY;
4366  }
4367 // std::cout<<"prob X1PY, X0PY, PXPY"<<probX1PY<<" "<<probX0PY<<" "<<probPXPY<<" "<<std::endl;
4368 // std::cout<<"gev,e0,e1,r0,r0: "<<gev<<" "<<e0<<" "<<e1<<" "<<r0<<" "<<r1<<std::endl;
4369 // std::cout<<"prob,pe0r0,pe0r1,pe1r0,pe1r1: "<<prob<<" "<<probe0r0<<" "<<probe0r1<<" "<<probe1r0<<" "<<probe1r1<<std::endl;
4370  if(interPos==false||probPXPY<0){
4371  unsigned int plane0 = int(int(plane)/std::abs(cos));
4372  unsigned int plane1 = int(int(plane)/std::abs(cos)) + 1;
4373  double r0 = double(plane/std::abs(cos))-double(plane0);
4374  double r1 = 1-r0;
4375  double probe0r0 = 1;
4376  double probe0r1 = 1;
4377  double probe1r0 = 1;
4378  double probe1r1 = 1;
4379  probe0r0 = RecoJMShower::GetPlaneDedxProb(shower, iReg, igev0, plane0, dedx, type);
4380  probe0r1 = RecoJMShower::GetPlaneDedxProb(shower, iReg, igev0, plane1, dedx, type);
4381  probe1r0 = RecoJMShower::GetPlaneDedxProb(shower, iReg, igev1, plane0, dedx, type);
4382  probe1r1 = RecoJMShower::GetPlaneDedxProb(shower, iReg, igev1, plane1, dedx, type);
4383  double probe0 = (r1*probe0r0+r0*probe0r1);
4384  double probe1 = (r1*probe1r0+r0*probe1r1);
4385  if(igev0!=igev1)probPXPY = (e1*probe0+e0*probe1)/(e0+e1);
4386  if(igev0==igev1)probPXPY = probe0;
4387  }
4388 
4389  return probPXPY;
4390  }
4391 /*
4392  double RecoJMShower::GetInterPlaneDedxProb(rb::Prong shower, unsigned int plane, double dedx, const int type, unsigned int nmip){
4393  int i = shower.ID();
4394  double gev = RecoJMShower::Energy(shower);
4395  double cos = shower.Dir()[2];
4396  double elower[11] = {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75};
4397  double eupper[11] = {0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25};
4398  double epoint[11];
4399  for(int ig=0;ig<11;ig++){
4400  epoint[ig] = (elower[ig]+eupper[ig])/2.;
4401  }
4402  int igev0 = 0;
4403  int igev1 = 0;
4404  if(gev>epoint[10]){igev0=10;igev1=10;}
4405  else if(gev<epoint[0]){igev0=0;igev1=0;}
4406  else{
4407  for(int ig=0;ig<10;ig++){
4408  if(gev>epoint[ig]&&gev<epoint[ig+1]){igev0=ig;igev1=ig+1;}
4409  }
4410  }
4411  if(igev0>10)igev0=10;
4412  if(igev1>10)igev1=10;
4413  unsigned int plane0 = int(int(plane)/std::abs(cos));
4414  unsigned int plane1 = int(int(plane)/std::abs(cos)) + 1;
4415  double r0 = double(plane/std::abs(cos))-double(plane0);
4416  double r1 = 1-r0;
4417  double e0 = gev-epoint[igev0];
4418  double e1 = epoint[igev1] - gev;
4419  double probe0r0 = 1;
4420  double probe0r1 = 1;
4421  double probe1r0 = 1;
4422  double probe1r1 = 1;
4423 
4424  if(plane1<=10){
4425  probe0r0 = RecoJMShower::GetPlaneDedxProb(shower, igev0, plane0, dedx, nmip, type);
4426  probe0r1 = RecoJMShower::GetPlaneDedxProb(shower, igev0, plane1, dedx, nmip, type);
4427  probe1r0 = RecoJMShower::GetPlaneDedxProb(shower, igev1, plane0, dedx, nmip, type);
4428  probe1r1 = RecoJMShower::GetPlaneDedxProb(shower, igev1, plane1, dedx, nmip, type);
4429 
4430 // if(int(jmshower::kElectron)==type&&shower.ID()==0&&fShowerPClusterCol[i].size()>20&&igev0>=3){
4431 // std::cout<<"e probe0r0, type, tot, mip "<<type<<" "<<totprobe0r0<<" "<<probe0r0<<std::endl;
4432 // std::cout<<"e probe0r1, type, tot, mip "<<type<<" "<<totprobe0r1<<" "<<probe0r1<<std::endl;
4433 // std::cout<<"e probe1r0, type, tot, mip "<<type<<" "<<totprobe1r0<<" "<<probe1r0<<std::endl;
4434 // std::cout<<"e probe1r1, type, tot, mip "<<type<<" "<<totprobe1r1<<" "<<probe1r1<<std::endl;
4435 // double probeE = RecoJMShower::GetPlaneDedxProb(shower, igev0, plane0, dedx, nmip, 0);
4436 // double probePi0 = RecoJMShower::GetPlaneDedxProb(shower, igev0, plane0, dedx, nmip, 3);
4437 // double totprobeE = RecoJMShower::GetPlaneDedxProb(shower, igev0, plane0, dedx, 0);
4438 // double totprobePi0 = RecoJMShower::GetPlaneDedxProb(shower, igev0, plane0, dedx, 3);
4439 // std::cout<<"================================= plane, igev0, nmip "<<plane0<<" "<<igev0<<" "<<nmip<<std::endl;
4440 // std::cout<<"dedx, mip e, pi0, delta === "<<dedx<<" "<<probeE<<" "<<probePi0<<" "<<probeE-probePi0<<std::endl;
4441 // std::cout<<"dedx, tot e, pi0, delta === "<<dedx<<" "<<totprobeE<<" "<<totprobePi0<<" "<<totprobeE-totprobePi0<<std::endl;
4442 // }
4443  }
4444  else{
4445  probe0r0 = RecoJMShower::GetPlaneDedxProb(shower, igev0, plane0, dedx, type);
4446  probe0r1 = RecoJMShower::GetPlaneDedxProb(shower, igev0, plane1, dedx, type);
4447  probe1r0 = RecoJMShower::GetPlaneDedxProb(shower, igev1, plane0, dedx, type);
4448  probe1r1 = RecoJMShower::GetPlaneDedxProb(shower, igev1, plane1, dedx, type);
4449  }
4450 
4451  double probe0 = (r1*probe0r0+r0*probe0r1);
4452  double probe1 = (r1*probe1r0+r0*probe1r1);
4453  double prob = 1;
4454  if(igev0!=igev1)prob = (e1*probe0+e0*probe1)/(e0+e1);
4455  if(igev0==igev1)prob = probe0;
4456 // std::cout<<"gev,e0,e1,r0,r0: "<<gev<<" "<<e0<<" "<<e1<<" "<<r0<<" "<<r1<<std::endl;
4457 // std::cout<<"prob,pe0r0,pe0r1,pe1r0,pe1r1: "<<prob<<" "<<probe0r0<<" "<<probe0r1<<" "<<probe1r0<<" "<<probe1r1<<std::endl;
4458  return prob;
4459  }
4460 */
4461 
4462  double RecoJMShower::GetDedxLongLL(rb::Prong shower, const int type){
4463  int i = shower.ID();
4464  TVector3 pos((shower.Start()[0]+RecoJMShower::GetShwStop(shower)[0])/2.0,(shower.Start()[1]+RecoJMShower::GetShwStop(shower)[1])/2.0,(shower.Start()[2]+RecoJMShower::GetShwStop(shower)[2])/2.0);
4465  int iReg = 0;
4466  if(pos[0]>=0&&pos[1]>=0)iReg = 0;
4467  if(pos[0]>=0&&pos[1]<0)iReg = 1;
4468  if(pos[0]<0&&pos[1]>=0)iReg = 2;
4469  if(pos[0]<=0&&pos[1]<=0)iReg = 3;
4470 
4471  double gev = RecoJMShower::Energy(shower);
4472  double elower[11] = {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75};
4473  double eupper[11] = {0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25};
4474  double epoint[11];
4475  for(int ig=0;ig<11;ig++){
4476  epoint[ig] = (elower[ig]+eupper[ig])/2.;
4477  }
4478  int igev0 = 0;
4479  int igev1 = 0;
4480  if(gev>epoint[10]){igev0=10;igev1=10;}
4481  else if(gev<epoint[0]){igev0=0;igev1=0;}
4482  else{
4483  for(int ig=0;ig<10;ig++){
4484  if(gev>epoint[ig]&&gev<epoint[ig+1]){igev0=ig;igev1=ig+1;}
4485  }
4486  }
4487  if(igev0>10)igev0=10;
4488  if(igev1>10)igev1=10;
4489  double e0 = gev-epoint[igev0];
4490  double e1 = epoint[igev1] - gev;
4491 
4492  double elowerm[41]= {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75,
4493  5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75,
4494  10.25,10.75,11.25,11.75,12.25,12.75,13.25,13.75,14.25,14.75,
4495  15.25,15.75,16.25,16.75,17.25,17.75,18.25,18.75,19.25,19.75
4496  };
4497 
4498  double eupperm[41]= { 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75,
4499  5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75,
4500  10.25,10.75,11.25,11.75,12.25,12.75,13.25,13.75,14.25,14.75,
4501  15.25,15.75,16.25,16.75,17.25,17.75,18.25,18.75,19.25,19.75, 20.25
4502  };
4503  double epointm[41];
4504  for(int ig=0;ig<41;ig++){
4505  epointm[ig] = (elowerm[ig]+eupperm[ig])/2.;
4506  }
4507 
4508  int igevm0 = 0;
4509  int igevm1 = 0;
4510  if(gev>epointm[40]){igevm0=40;igevm1=40;}
4511  else if(gev<epointm[0]){igevm0=0;igevm1=0;}
4512  else{
4513  for(int ig=0;ig<40;ig++){
4514  if(gev>epointm[ig]&&gev<epointm[ig+1]){igevm0=ig;igevm1=ig+1;}
4515  }
4516  }
4517  if(igevm0>40)igevm0=40;
4518  if(igevm1>40)igevm1=40;
4519  double em0 = gev-epointm[igevm0];
4520  double em1 = epointm[igevm1] - gev;
4521 
4522 
4523  double sumlgprob = 0;
4524 // unsigned int totplane = 90;
4525  unsigned int totplane = fShowerPClusterCol[i].size();
4526  unsigned int nonzeroplane = 0;
4527 
4528 /*
4529  unsigned int shnp = fShowerPClusterCol[i].size();
4530  if(int(jmshower::kElectron)==type)totplane=std::max(shnp,fHtExpPlaneE[iReg][igev1]);
4531  if(int(jmshower::kGamma)==type)totplane=std::max(shnp,fHtExpPlaneG[iReg][igev1]);
4532  if(int(jmshower::kMuon)==type)totplane=std::max(shnp,fHtExpPlaneMu[iReg][igev1]);
4533  if(int(jmshower::kPi0)==type)totplane=std::max(shnp,fHtExpPlanePi0[iReg][igev1]);
4534  if(int(jmshower::kHad)==type)totplane=std::max(shnp,fHtExpPlaneHad[iReg][igev1]);
4535  if(int(jmshower::kProton)==type)totplane=std::max(shnp,fHtExpPlaneP[iReg][igev1]);
4536  if(int(jmshower::kNeutron)==type)totplane=std::max(shnp,fHtExpPlaneN[iReg][igev1]);
4537  if(int(jmshower::kPion)==type)totplane=std::max(shnp,fHtExpPlanePi[iReg][igev1]);
4538  if(int(jmshower::kElectronQE)==type)totplane=std::max(shnp,fHtExpPlaneEqe[iReg][igev1]);
4539  if(int(jmshower::kElectronRES)==type)totplane=std::max(shnp,fHtExpPlaneEres[iReg][igev1]);
4540  if(int(jmshower::kElectronDIS)==type)totplane=std::max(shnp,fHtExpPlaneEdis[iReg][igev1]);
4541  if(int(jmshower::kElectronCOH)==type)totplane=std::max(shnp,fHtExpPlaneEcoh[iReg][igev1]);
4542 */
4543 
4544 
4545  //std::cout<<"type, energy, cos, np, totplane "<<type<<" "<<gev<<" "<<shower.Dir()[2]<<" "<<shnp<<" "<<totplane<<std::endl;
4546 
4547  for(unsigned int plane=0;plane<totplane;plane++){
4548 // TVector3 disttoedge = RecoJMShower::GetTrkPlaneDistToEdge(fShowerPClusterCol[i][plane].Plane(), shower);
4549 // double planerad = RecoJMShower::Radius(fShowerPClusterCol[i][plane], shower);
4550 // if(disttoedge.x()<2*planerad)continue;
4551 // if(disttoedge.y()<2*planerad)continue;
4552 // if(disttoedge.z()<2*planerad)continue;
4553 // double path = RecoJMShower::GetTrkHitPath(fShowerPClusterCol[i][plane].Plane(), 0, shower);
4554 // double ep = RecoJMShower::Energy(fShowerPClusterCol[i][plane],shower.ID());
4555 // double dedx = ep/path;
4556  double dedx = 0;
4557  if(plane<fShowerPClusterCol[i].size()){
4558  dedx = RecoJMShower::GetPlaneDedx(shower, plane);
4559  }
4560  double prob = 0;
4561 // unsigned int nmip = RecoJMShower::NMIPPlane(shower);
4562  if(type != int(jmshower::kElectronSG)){
4563  prob = RecoJMShower::GetInterPlaneDedxProb(shower, plane, dedx, type, pos, iReg, igev0, igev1, e0, e1);
4564  }else{
4565  prob = RecoJMShower::GetInterPlaneDedxProb(shower, plane, dedx, type, pos, iReg, igevm0, igevm1, em0, em1);
4566  }
4567 // prob = RecoJMShower::GetPlaneDedxProb(shower, i, hiteweightmap, 2,plane, dedx, type);
4568  if(prob<0.0001)prob=0.0001;
4569  double lgprob = log(prob);
4570  if(pBadChannel==true){ //Only consider non-zero plane
4571  if(dedx>0.0000001){
4572  nonzeroplane++;
4573  sumlgprob+=lgprob;
4574  }
4575  }else{
4576  nonzeroplane++;
4577  sumlgprob+=lgprob;
4578  }
4579  //std::cout<<"type: "<<type<<" "<<plane<<" "<<prob<<" "<<sumlgprob<<std::endl;
4580  }
4581  if(pBadChannel==true){
4582  if(nonzeroplane>0)sumlgprob = sumlgprob/(double)nonzeroplane;
4583  }else{
4584  if(totplane>0)sumlgprob = sumlgprob/(double)totplane;
4585  }
4586 // std::cout<<"sumlgprob, totplane "<<sumlgprob<<" "<<totplane<<std::endl;
4587 // sumlgprob = sumlgprob/double(fShowerPClusterCol[i].size());
4588  if(sumlgprob>999)sumlgprob=999;
4589  if(sumlgprob<-999)sumlgprob=-999;
4590  //std::cout<<"DedxLongLL type: "<<type<<" ll "<<sumlgprob<<std::endl;
4591  return sumlgprob;
4592  }
4593 
4594  double RecoJMShower::GetDedxInvLongLL(rb::Prong shower, const int type){
4595  int i = shower.ID();
4596  TVector3 pos((shower.Start()[0]+RecoJMShower::GetShwStop(shower)[0])/2.0,(shower.Start()[1]+RecoJMShower::GetShwStop(shower)[1])/2.0,(shower.Start()[2]+RecoJMShower::GetShwStop(shower)[2])/2.0);
4597  int iReg = 0;
4598  if(pos[0]>=0&&pos[1]>=0)iReg = 0;
4599  if(pos[0]>=0&&pos[1]<0)iReg = 1;
4600  if(pos[0]<0&&pos[1]>=0)iReg = 2;
4601  if(pos[0]<=0&&pos[1]<=0)iReg = 3;
4602 
4603  double gev = RecoJMShower::Energy(shower);
4604  double elower[11] = {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75};
4605  double eupper[11] = {0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25};
4606  double epoint[11];
4607  for(int ig=0;ig<11;ig++){
4608  epoint[ig] = (elower[ig]+eupper[ig])/2.;
4609  }
4610  int igev0 = 0;
4611  int igev1 = 0;
4612  if(gev>epoint[10]){igev0=10;igev1=10;}
4613  else if(gev<epoint[0]){igev0=0;igev1=0;}
4614  else{
4615  for(int ig=0;ig<10;ig++){
4616  if(gev>epoint[ig]&&gev<epoint[ig+1]){igev0=ig;igev1=ig+1;}
4617  }
4618  }
4619  if(igev0>10)igev0=10;
4620  if(igev1>10)igev1=10;
4621  double e0 = gev-epoint[igev0];
4622  double e1 = epoint[igev1] - gev;
4623 
4624 
4625  double sumlgprob = 0;
4626 // unsigned int totplane = 90;
4627  unsigned int totplane = fShowerPClusterCol[i].size();
4628  unsigned int nonzeroplane = 0;
4629 
4630 /*
4631  unsigned int shnp = fShowerPClusterCol[i].size();
4632  if(int(jmshower::kElectron)==type)totplane=std::max(shnp,fHtExpPlaneE[iReg][igev1]);
4633  if(int(jmshower::kGamma)==type)totplane=std::max(shnp,fHtExpPlaneG[iReg][igev1]);
4634  if(int(jmshower::kMuon)==type)totplane=std::max(shnp,fHtExpPlaneMu[iReg][igev1]);
4635  if(int(jmshower::kPi0)==type)totplane=std::max(shnp,fHtExpPlanePi0[iReg][igev1]);
4636  if(int(jmshower::kHad)==type)totplane=std::max(shnp,fHtExpPlaneHad[iReg][igev1]);
4637  if(int(jmshower::kProton)==type)totplane=std::max(shnp,fHtExpPlaneP[iReg][igev1]);
4638  if(int(jmshower::kNeutron)==type)totplane=std::max(shnp,fHtExpPlaneN[iReg][igev1]);
4639  if(int(jmshower::kPion)==type)totplane=std::max(shnp,fHtExpPlanePi[iReg][igev1]);
4640  if(int(jmshower::kElectronQE)==type)totplane=std::max(shnp,fHtExpPlaneEqe[iReg][igev1]);
4641  if(int(jmshower::kElectronRES)==type)totplane=std::max(shnp,fHtExpPlaneEres[iReg][igev1]);
4642  if(int(jmshower::kElectronDIS)==type)totplane=std::max(shnp,fHtExpPlaneEdis[iReg][igev1]);
4643  if(int(jmshower::kElectronCOH)==type)totplane=std::max(shnp,fHtExpPlaneEcoh[iReg][igev1]);
4644 */
4645 
4646 
4647  for(unsigned int plane=0;plane<totplane;plane++){
4648 // TVector3 disttoedge = RecoJMShower::GetTrkPlaneDistToEdge(fShowerPClusterCol[i][plane].Plane(), shower);
4649 // double planerad = RecoJMShower::Radius(fShowerPClusterCol[i][plane], shower);
4650 // if(disttoedge.x()<2*planerad)continue;
4651 // if(disttoedge.y()<2*planerad)continue;
4652 // if(disttoedge.z()<2*planerad)continue;
4653 // double path = RecoJMShower::GetTrkHitPath(fShowerPClusterCol[i][plane].Plane(), 0, shower);
4654 // double ep = RecoJMShower::Energy(fShowerPClusterCol[i][plane],shower.ID());
4655 // double dedx = ep/path;
4656  double dedx = 0;
4657  if(totplane-plane<fShowerPClusterCol[i].size()){
4658  dedx = RecoJMShower::GetPlaneDedx(shower, totplane-plane);
4659  }
4660  double prob = 0;
4661 // unsigned int nmip = RecoJMShower::NMIPPlane(shower);
4662  prob = RecoJMShower::GetInterPlaneDedxProb(shower, plane, dedx, type, pos, iReg, igev0, igev1, e0, e1);
4663 // prob = RecoJMShower::GetPlaneDedxProb(shower, i, hiteweightmap, 2,plane, dedx, type);
4664  if(prob<0.0001)prob=0.0001;
4665  double lgprob = log(prob);
4666  if(pBadChannel==true){ //Only consider non-zero plane
4667  if(dedx>0.0000001){
4668  nonzeroplane++;
4669  sumlgprob+=lgprob;
4670  }
4671  }else{
4672  nonzeroplane++;
4673  sumlgprob+=lgprob;
4674  }
4675  }
4676  if(pBadChannel==true){
4677  if(nonzeroplane>0)sumlgprob = sumlgprob/(double)nonzeroplane;
4678  }else{
4679  if(totplane>0)sumlgprob = sumlgprob/(double)totplane;
4680  }
4681 // sumlgprob = sumlgprob/double(fShowerPClusterCol[i].size());
4682  if(sumlgprob>999)sumlgprob=999;
4683  if(sumlgprob<-999)sumlgprob=-999;
4684  return sumlgprob;
4685  }
4686 
4687 /*
4688  double RecoJMShower::GetDedxLongLL(rb::Prong shower, const int type, unsigned int startp, int contp, art::Event const& evt){
4689  int i = shower.ID();
4690  double sumlgprob = 0;
4691  if(int(fShowerPClusterCol[i].size()-startp)<contp)return -999;
4692  unsigned int newstartp = RecoJMShower::ContStartPlane(shower, startp, contp);
4693 
4694  for(unsigned int plane= newstartp; plane<fShowerPClusterCol[i].size();plane++){
4695  TVector3 disttoedge = RecoJMShower::GetTrkPlaneDistToEdge(fShowerPClusterCol[i][plane].Plane(), shower);
4696 // double planerad = RecoJMShower::Radius(fShowerPClusterCol[i][plane], shower);
4697 // if(disttoedge.x()<2*planerad)continue;
4698 // if(disttoedge.y()<2*planerad)continue;
4699 // if(disttoedge.z()<2*planerad)continue;
4700 // double path = RecoJMShower::GetTrkHitPath(fShowerPClusterCol[i][plane].Plane(), 0, shower);
4701 // double ep = RecoJMShower::Energy(fShowerPClusterCol[i][plane],shower.ID());
4702 // double dedx = ep/path;
4703  double dedx = RecoJMShower::GetPlaneDedx(shower, plane);
4704  double prob = 0;
4705 // unsigned int nmip = RecoJMShower::NMIPPlane(shower, newstartp);
4706  if(newstartp==startp)prob = RecoJMShower::GetInterPlaneDedxProb(shower, plane, dedx, type);
4707  if(newstartp>startp)prob = RecoJMShower::GetInterPlaneDedxProb(shower, plane-newstartp, dedx, type);
4708 // prob = RecoJMShower::GetPlaneDedxProb(shower, i, hiteweightmap, 2,plane, dedx, type);
4709  if(prob<0.0001)prob = 0.0001;
4710  double lgprob = log(prob);
4711  sumlgprob+=lgprob;
4712 // std::cout<<"i, plane, dedx, prob "<<i<<" "<<plane<<" "<<dedx<<" "<<prob<<std::endl;
4713  }
4714  sumlgprob = sumlgprob/double(fShowerPClusterCol[i].size() - newstartp);
4715  return sumlgprob;
4716  if(sumlgprob>999)sumlgprob=999;
4717  if(sumlgprob<-999)sumlgprob=-999;
4718  }
4719 */
4720 
4721 /*
4722  double RecoJMShower::GetCellTransDedx(rb::Prong shower, unsigned int cell){
4723  int i = shower.ID();
4724  double totpath = 0 ;
4725  for(unsigned int plane=0;plane<fShowerPClusterCol[i].size();plane++){
4726  double path = RecoJMShower::GetTrkHitPath(fShowerPClusterCol[i][plane].Plane(), 0, shower);
4727  totpath+=path;
4728  }
4729  double ep = RecoJMShower::DepositEnergy(fShowerTCClusterCol[i][cell], i);
4730  double dedx = ep/totpath;
4731  if(isrealdata==true&&pCorDataDedx)dedx=dedx/RecoJMShower::CellTransDataMC(dedx, fDetid);
4732  return dedx;
4733  }
4734 */
4735 
4736  double RecoJMShower::GetCellTransDedx(rb::Prong shower, unsigned int cell){// calculate transverse cell dE/dx based on centroid
4737 //if(cell!=0)return 0;
4738  int i = shower.ID();
4739  double totpath = 0 ;
4740  double ep = 0;
4741  for(unsigned int plane=0;plane<fShowerPClusterCol[i].size();plane++){
4742  unsigned int absplane = fShowerPClusterCol[i][plane].Plane();
4743  art::ServiceHandle<geo::Geometry> geom;
4744  double sumPos = 0;
4745  double gev = 0;
4746  for (unsigned int nh=0;nh<fShowerPClusterCol[i][plane].NCell();nh++) {
4747  art::Ptr<rb::CellHit> cellhit = fShowerPClusterCol[i][plane].Cell(nh);
4748  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
4749  double cellEnergy = fHitEWeightMap[key][i];
4750  double cellPos = (double)cellhit->Cell() + 0.5;
4751  sumPos += cellEnergy*cellPos;
4752  gev += cellEnergy;
4753  }
4754  if(gev<0.00001)continue;
4755  double avgPos = sumPos/gev;
4756  int cid = (int)avgPos;// calculate centroid
4757  double wm = avgPos - (double)cid;//energy splitting for the centroid cell
4758  double wp = 1-wm; //energy splitting for the centroid cell
4759  if(cid<0)cid=0;
4760  if(cid>(int)geom->Plane(absplane)->Ncells()-1)cid=(int)geom->Plane(absplane)->Ncells()-1;
4761  double dir = 1.;
4762  if(std::abs(shower.Dir()[2])>0.000001)dir=std::abs(shower.Dir()[2]);
4763 // double dx=shower.Dir()[0];// Validate dir
4764 // double dy=shower.Dir()[1];
4765 // double dz=shower.Dir()[2];
4766 // if(geom->Plane(absplane)->View()==geo::kX&&dx*dx+dz*dz>1e-10)dir=std::abs(dz/sqrt(dx*dx+dz*dz));
4767 // if(geom->Plane(absplane)->View()==geo::kY&&dy*dy+dz*dz>1e-10)dir=std::abs(dz/sqrt(dy*dy+dz*dz));
4768  int ic = (int)cell + 1;
4769  double dph = ((double)ic/dir)+wm;//upper limit for cells included in a transeverse cell index
4770  double dpl = ((double)(ic-1)/dir)+wm;//lower limit for cells included in a transeverse cell index
4771  if(cid+(int)dph<(int)geom->Plane(plane)->Ncells()){
4772  double eps1 = dph-(int)dph;
4773  double eps2 = (int)dpl+1-dpl;
4774  geo::OfflineChan keye1(absplane, cid+(int)dph);
4775  geo::OfflineChan keye2(absplane, cid+(int)dpl);
4776  double e1 = 0;
4777  double e2 = 0;
4778  double e = 0;
4779  if((int)fHitEWeightMap[keye1].size()>i)e1 = eps1*fHitEWeightMap[keye1][i];
4780  if((int)fHitEWeightMap[keye2].size()>i)e2 = eps2*fHitEWeightMap[keye2][