ParticleIDAlg.cxx
Go to the documentation of this file.
1 
2 ////////////////////////////////////////////////////////////////////////
3 /// \file ParticleIDAlg.cxx
4 /// \brief Calculate dE/dx and log likelihoods for different particle hypotheses.
5 /// This information will be of use for identification of any particle that
6 /// uses rb::Showers as a starting point.
7 ///
8 /// \todo Remove hardcoded physical constants.
9 ///
10 ///
11 /// \version $Id: ParticleIDAlg.cxx,v 1.11 2012-10-26 22:19:03 bckhouse Exp $
12 /// \author Alex Smith (smith@physics.umn.edu)
13 ////////////////////////////////////////////////////////////////////////
14 
15 
17 #include "TLorentzVector.h"
20 #include "GeometryObjects/Geo.h"
21 
22 #include <iterator>
23 #include <algorithm>
24 #include <map>
25 #include <cmath>
26 
27 namespace slid {
28 
29  //*********************************************************************************************
32  NuEEnergyAlg* NuEEnergyAlgIn,
33  bool loadLibs)
34  : fClosestPi0Mass(-9999),
35  fDedxDistributions(slid::kDEDXPARTICLETYPESIZE),
36  fShower(0),
37  fNuEEnergyAlg(NuEEnergyAlgIn),
38  fAsym(1.0),
39  fMipRangeHigh(0.005),
40  fMipRangeLow(0.0005),
41  fDetId(detId)
42  {
43  if (loadLibs){
44  LoadDedxHistogramFiles(pset, detId);
45  }
46  }
47 
48  //*********************************************************************************************
50  {
51  }
52 
53  //*********************************************************************************************
56  {
57  }
58 
59  //*********************************************************************************************
60  TFile * LoadHistFromFile(TString filePath)
61  {
62  std::cout << "LID Dedx File " << filePath << " is loaded." << std::endl;
63  return new TFile(filePath);
64  }
65 
66  /// \brief Calculate transverse dE/dx log-likelihood for specific particle hypothesis
67  //*********************************************************************************************
69  const rb::Shower* vShower,
70  const std::vector< const rb::Shower* > showercol,
71  const art::Event& evt)
72  {
73  // Lazy setting of shower
74  SetShower(vShower, showercol, evt);
75  std::map< int, float>::iterator llIter = fPartTransLL.find(partHyp);
76  if (llIter == fPartTransLL.end()) {
77  double sumLogProb = 0;
78  int transPlane = 0;
79  for (int iTransPlane = 0; iTransPlane < kNumTransversePlane; iTransPlane++) {
80  double dedx = CalcPlaneTransverseDedx(iTransPlane);
81  if (dedx < -0.5) continue;
82  double prob = 0;
83  prob = CalcInterCellTransDedxProb(partHyp, iTransPlane);
84  if (prob < 0.0001) prob = 0.0001;
85  double logProb = log(prob);
86  sumLogProb += logProb;
87  transPlane++;
88  }
89  if (transPlane >0 ) sumLogProb = sumLogProb / (double)transPlane;
90  if (sumLogProb > 999)sumLogProb = 999;
91  if (sumLogProb<-999)sumLogProb = -999;
92 
93  fPartTransLL[partHyp] = sumLogProb;
94  return sumLogProb;
95  }
96  return llIter->second;
97  }
98 
99  /// \brief Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis
100  //*********************************************************************************************
102  const rb::Shower* vShower,
103  const std::vector< const rb::Shower* > showercol,
104  const art::Event& evt)
105  {
106  // Lazy setting of shower
107  SetShower(vShower, showercol, evt);
108 
109  std::map< int, float>::iterator llIter = fPartLongLL.find(partHyp);
110  if (llIter == fPartLongLL.end()) {
111  double sumLogLikelihood = 0;
112  unsigned int totplane = fPlaneHits.size();
113  //unsigned int nonzeroplane = 0;
114  unsigned int planecount = 0;
115 
116  if (vShower->Dir()[2] >= 0.0){ //forward going shower
117  for (std::map<int, std::map<int, int> >::iterator itr = fPlaneHits.begin(); itr != fPlaneHits.end(); ++itr){
118  double prob = 0;
119  prob = CalcInterPlaneDedxProb(partHyp, planecount, itr->first);
120  if (prob < 0.0001)prob = 0.0001;
121  double lgprob = log(prob);
122  //nonzeroplane++;
123  sumLogLikelihood += lgprob;
124  planecount++;
125  }
126  }
127  else{ //backward going shower
128  for (std::map<int, std::map<int, int> >::reverse_iterator itr = fPlaneHits.rbegin(); itr != fPlaneHits.rend(); ++itr){
129  double prob = 0;
130  prob = CalcInterPlaneDedxProb(partHyp, planecount, itr->first);
131 
132  if (prob < 0.0001)prob = 0.0001;
133  double lgprob = log(prob);
134  //nonzeroplane++;
135  sumLogLikelihood += lgprob;
136  planecount++;
137  }
138  }
139  if (totplane > 0) sumLogLikelihood = sumLogLikelihood / (double) totplane;
140  if (sumLogLikelihood > 999) sumLogLikelihood = 999;
141  if (sumLogLikelihood<-999) sumLogLikelihood = -999;
142  fPartLongLL[partHyp] = sumLogLikelihood;
143  return sumLogLikelihood;
144  }
145  return llIter->second;
146  }
147 
148  /// \brief Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis
149  //*********************************************************************************************
151  const rb::Shower* vShower,
152  const std::vector< const rb::Shower* > showercol,
153  const art::Event& evt)
154  {
155  rb::Shower newShower = *vShower;
156  newShower.SetStart( vShower->Stop());
157  newShower.SetDir( -(vShower->Dir()) );
158  return DedxLongLL(partHyp, &(newShower), showercol, evt);
159  }
160 
161  /**
162  * \brief Set rb::Shower object to be analyzed. Calculate its total energy, which will be needed
163  * for calculations.
164  *
165  * @param vShower
166  * @param evt
167  */
168  //*********************************************************************************************
169  void ParticleIDAlg::SetShower(const rb::Shower* vShower,
170  const std::vector< const rb::Shower* > showercol,
171  const art::Event& evt)
172  {
173  //
174  // Only recalculate shower quantities if a different event or shower is presented or if
175  // fShower is not already set
176  //
177 
178  if ( fShower != NULL &&
179  evt.id() == eventID &&
180  fShowerCol.size() == showercol.size() &&
181  std::equal( fShowerCol.begin(), fShowerCol.end(), showercol.begin()) &&
182  fShower->ID() == vShower->ID())
183  return;
184 
185  // Save to compare next one
186  fShower = vShower;
187  fShowerCol = showercol;
188  eventID = evt.id();
189 
190  //clear map of plane index to global plane
191  fMapIdxToGlobalPlane.clear();
192 
193  // Clear cache of longitudinal dedx values by longitudinal and tranverse plane
194  fDedxLongByPlane.clear();
195  fDedxTransByPlane.clear();
196 
197  fPartLongLL.clear();
198  fPartTransLL.clear();
199 
200  fTransEnergy.clear();
201  fTransPathLength.clear();
202  fPi0Index = -5;
203 
204  fIxyz.SetXYZ(0.0,0.0,0.0);
205  fAsym = 1.0;
206 
207  // Calculate and save the total energy of the shower in a member field
210 
211 
212  // Fill a map where the first element is
213  // plane number and the second is vector of global
214  // cell indices in a shower of cells that belong to
215  // that plane.
216  fPlaneHits = CalcPlaneHits( vShower);
217 
218  // Fill a map where the first index is the
219  // number of planes away from the track core, and the
220  // second is a vector of global cell indices in a
221  // shower that are that many planes from the track
222  // core. It collapses the transverse profile of the
223  // shower to a single plane.
224  fPlaneCentroid = CalcPlaneCentroid( vShower );
225 
226  //
227  // Determine detector XY region index
228  //
230  //
231  // Determine energy bins
232  //
234 
235  //calculate shower asymmetry and ineria
236  CalcAsymIneria();
237 
238  }
239 
240  //*********************************************************************************************
241  std::map< int, std::map< int, int> > ParticleIDAlg::CalcPlaneHits(const rb::Shower* vShower)
242  {
243  std::map< int, std::map< int, int> > mapP;
244 
245  std::map<int, int> dummymap;
246 
247  int startPlane = vShower->MinPlane();
248  int stopPlane = vShower->MaxPlane();
249  for (int i=startPlane; i<=stopPlane; ++i){
250  mapP[i] = dummymap;
251  }
252 
253  int ncells = vShower->NCell();
254  for (int i = 0; i < ncells; ++i) {
255 
256  mapP[ vShower->Cell(i)->Plane() ][vShower->Cell(i)->Cell()] = i;
257  }// end loop over cells
258 
259  //now map the plane index with the global plane number
260  unsigned int planecount = 0;
261  if (vShower->Dir()[2] >= 0.0){ //forward going shower
262  for (std::map<int, std::map<int, int> >::iterator itr = mapP.begin(); itr != mapP.end(); ++itr){
263  fMapIdxToGlobalPlane[planecount] = itr->first;
264  planecount++;
265  }
266  }
267  else{ //backward going shower
268  for (std::map<int, std::map<int, int> >::reverse_iterator itr = mapP.rbegin(); itr != mapP.rend(); ++itr){
269  fMapIdxToGlobalPlane[planecount] = itr->first;
270  planecount++;
271  }
272  }
273 
274  return mapP;
275  }// end of CalcPlaneHits
276 
277  //*********************************************************************************************
278  std::map< int, double > ParticleIDAlg::CalcPlaneCentroid(const rb::Shower* vShower)
279  {
280 
281  std::map< int, double > mapT;
282 
283  //loop plane by plane to calculate weighted cell center and then transverse cells
284  for(std::map<int, std::map<int, int> >::iterator itr=fPlaneHits.begin(); itr!=fPlaneHits.end(); ++itr){
285  std::map<int, int> planecells = itr->second;
286  //loop over hits in plane to determine energy weighted center
287  double sumPos = 0.0;
288  double planeGev = 0.0;
289  for(std::map<int, int>::iterator itrC = planecells.begin(); itrC!=planecells.end(); ++itrC){
290  double cellPos = itrC->first + 0.5;
291  double cellE = fNuEEnergyAlg->CellEnergy(vShower, fShowerCol, itrC->second, eventID );
292  if ((cellE == -5) || std::isnan(cellE)) cellE = 0.0;
293  sumPos += cellPos*cellE;
294  planeGev += cellE;
295  }
296  double avgPos = -999.0; //non physical cell position as the default
297  if (planeGev > 0.00001) avgPos = sumPos/planeGev;
298  mapT[itr->first] = avgPos;
299  }
300 
301  return mapT;
302  }
303 
304 
305  //*********************************************************************************************
307  novadaq::cnv::DetId detectorID)
308  {
309 
310  std::string libPath = util::EnvExpansion(pset.get< std::string >("LibraryPath"));
311  std::string thistname = pset.get< std::string >("HistFDDedxElectron");
312  std::cout << "Loading LID library file from " << libPath << std::endl;
313  if (detectorID == novadaq::cnv::DetId::kFARDET) {
314 
315  /// FD dE/dx input histogram files
316  std::cout<<"---- dE/dx histogram file for electron\n";
317  fDedxDistributions[slid::kELECTRON].Initialize(pset, "HistFDDedxElectron");
318  std::cout<<"---- dE/dx histogram file for photon\n";
319  fDedxDistributions[slid::kPHOTON].Initialize(pset, "HistFDDedxPhoton");
320  std::cout<<"---- dE/dx histogram file for muon\n";
321  fDedxDistributions[slid::kMUON].Initialize(pset, "HistFDDedxMuon");
322  std::cout<<"---- dE/dx histogram file for pi0\n";
323  fDedxDistributions[slid::kPI0].Initialize(pset, "HistFDDedxPi0");
324  std::cout<<"---- dE/dx histogram file for proton\n";
325  fDedxDistributions[slid::kPROTON].Initialize(pset, "HistFDDedxProton");
326  std::cout<<"---- dE/dx histogram file for neutron\n";
327  fDedxDistributions[slid::kNEUTRON].Initialize(pset, "HistFDDedxNeutron");
328  std::cout<<"---- dE/dx histogram file for charged pion\n";
329  fDedxDistributions[slid::kPION].Initialize(pset, "HistFDDedxPion");
330  std::cout<<"---- dE/dx histogram file for nue CCQE\n";
331  fDedxDistributions[slid::kELECTRONCCQE].Initialize(pset, "HistFDDedxElectronCCQE");
332  std::cout<<"---- dE/dx histogram file for nue CCRes\n";
333  fDedxDistributions[slid::kELECTRONCCRES].Initialize(pset, "HistFDDedxElectronCCRES");
334  std::cout<<"---- dE/dx histogram file for nue CC DIS\n";
335  fDedxDistributions[slid::kELECTRONCCDIS].Initialize(pset, "HistFDDedxElectronCCDIS");
336  std::cout<<"---- dE/dx histogram file for nue Coherent\n";
337  fDedxDistributions[slid::kELECTRONCCCOH].Initialize(pset, "HistFDDedxElectronCCCOH");
338 
339  } else if (detectorID == novadaq::cnv::DetId::kNDOS) {
340 
341  /// NDOS dE/dx input histograms
342 
343  std::cout<<"---- dE/dx histogram file for electron\n";
344  fDedxDistributions[slid::kELECTRON].Initialize(pset, "HistNDOSDedxElectron");
345  std::cout<<"---- dE/dx histogram file for photon\n";
346  fDedxDistributions[slid::kPHOTON].Initialize(pset, "HistNDOSDedxPhoton");
347  std::cout<<"---- dE/dx histogram file for muon\n";
348  fDedxDistributions[slid::kMUON].Initialize(pset, "HistNDOSDedxMuon");
349  std::cout<<"---- dE/dx histogram file for pi0\n";
350  fDedxDistributions[slid::kPI0].Initialize(pset, "HistNDOSDedxPi0");
351  std::cout<<"---- dE/dx histogram file for proton\n";
352  fDedxDistributions[slid::kPROTON].Initialize(pset, "HistNDOSDedxProton");
353  std::cout<<"---- dE/dx histogram file for neutron\n";
354  fDedxDistributions[slid::kNEUTRON].Initialize(pset, "HistNDOSDedxNeutron");
355  std::cout<<"---- dE/dx histogram file for charged pion\n";
356  fDedxDistributions[slid::kPION].Initialize(pset, "HistNDOSDedxPion");
357  std::cout<<"---- dE/dx histogram file for nue CCQE\n";
358  fDedxDistributions[slid::kELECTRONCCQE].Initialize(pset, "HistNDOSDedxElectronCCQE");
359  std::cout<<"---- dE/dx histogram file for nue CCRes\n";
360  fDedxDistributions[slid::kELECTRONCCRES].Initialize(pset, "HistNDOSDedxElectronCCRES");
361  std::cout<<"---- dE/dx histogram file for nue CC DIS\n";
362  fDedxDistributions[slid::kELECTRONCCDIS].Initialize(pset, "HistNDOSDedxElectronCCDIS");
363  std::cout<<"---- dE/dx histogram file for nue Coherent\n";
364  fDedxDistributions[slid::kELECTRONCCCOH].Initialize(pset, "HistNDOSDedxElectronCCCOH");
365 
366  } else if (detectorID == novadaq::cnv::DetId::kNEARDET) {
367 
368  /// ND dE/dx input histograms
369 
370  std::cout<<"---- dE/dx histogram file for electron\n";
371  fDedxDistributions[slid::kELECTRON].Initialize(pset, "HistNDDedxElectron");
372  std::cout<<"---- dE/dx histogram file for photon\n";
373  fDedxDistributions[slid::kPHOTON].Initialize(pset, "HistNDDedxPhoton");
374  std::cout<<"---- dE/dx histogram file for muon\n";
375  fDedxDistributions[slid::kMUON].Initialize(pset, "HistNDDedxMuon");
376  std::cout<<"---- dE/dx histogram file for pi0\n";
377  fDedxDistributions[slid::kPI0].Initialize(pset, "HistNDDedxPi0");
378  std::cout<<"---- dE/dx histogram file for proton\n";
379  fDedxDistributions[slid::kPROTON].Initialize(pset, "HistNDDedxProton");
380  std::cout<<"---- dE/dx histogram file for neutron\n";
381  fDedxDistributions[slid::kNEUTRON].Initialize(pset, "HistNDDedxNeutron");
382  std::cout<<"---- dE/dx histogram file for charged pion\n";
383  fDedxDistributions[slid::kPION].Initialize(pset, "HistNDDedxPion");
384  std::cout<<"---- dE/dx histogram file for nue CCQE\n";
385  fDedxDistributions[slid::kELECTRONCCQE].Initialize(pset, "HistNDDedxElectronCCQE");
386  std::cout<<"---- dE/dx histogram file for nue CC Res\n";
387  fDedxDistributions[slid::kELECTRONCCRES].Initialize(pset, "HistNDDedxElectronCCRES");
388  std::cout<<"---- dE/dx histogram file for nue CC DIS\n";
389  fDedxDistributions[slid::kELECTRONCCDIS].Initialize(pset, "HistNDDedxElectronCCDIS");
390  std::cout<<"---- dE/dx histogram file for nue Coherent\n";
391  fDedxDistributions[slid::kELECTRONCCCOH].Initialize(pset, "HistNDDedxElectronCCCOH");
392  } else {
393  throw cet::exception("DetectorTypeNotDefined")
394  << "Detector type " << detectorID
395  << " not recognized. See LIDBuilder_module.cc for list of defined algorithms. "
396  << "Only novadaq::cnv::DetId::kNEARDET, kFARDET, kNDOS defined.";
397 
398  }
399  return true;
400  }
401 
402  // bool ParticleIDAlg::LoadDedxHistograms() {
403  //
404  // for (int iParticleType = 0; iParticleType < slid::kDEDXPARTICLETYPESIZE;
405  // iParticleType++) {
406  // fDedxDistributions[iParticleType].LoadDedxHistograms();
407  // }
408  // return true;
409  // }
410 
411  ///
412  /// \brief Determine detector XY region
413  ///
414 
415  //*********************************************************************************************
417  {
418 
419  TVector3 pos((fShower->Start()[0] + fShower->Stop()[0]) / 2.0,
420  (fShower->Start()[1] + fShower->Stop()[1]) / 2.0,
421  (fShower->Start()[2] + fShower->Stop()[2]) / 2.0);
422  iReg = 0;
423  //only switch quadrants if this is the Far Detector, for ND whole detector is ~quadrant 0
424  // (near the readout in both views)
426  if (pos[0] >= 0 && pos[1] >= 0)iReg = 0;
427  if (pos[0] >= 0 && pos[1] < 0)iReg = 1;
428  if (pos[0] < 0 && pos[1] >= 0)iReg = 2;
429  if (pos[0] <= 0 && pos[1] <= 0)iReg = 3;
430  }
431  return;
432  }
433 
434  ///
435  /// \brief Calculate energy bin indices and energy values from the histograms
436  /// corresponding energy of shower to be used to interpolate between bins in energy
437  ///
438 
439  //*********************************************************************************************
441  {
442  // int& igev0, int& igev1,
443  // double& e0, double& e1
444 
445  double energyBinBounds[slid::kNumEnergyBin];
446  for (int ig = 0; ig < slid::kNumEnergyBin; ig++) {
447  energyBinBounds[ig] = (slid::kEnergyBins[ig] + kEnergyBins[ig + 1]) / 2.;
448  }
449 
450  igev0 = 0;
451  igev1 = 0;
452  if (fShowerEnergyGev > energyBinBounds[slid::kNumEnergyBin - 1]) {
453  igev0 = slid::kNumEnergyBin - 1;
454  igev1 = slid::kNumEnergyBin - 1;
455  } else if (fShowerEnergyGev < energyBinBounds[0]) {
456  igev0 = 0;
457  igev1 = 0;
458  } else {
459  for (int ig = 0; ig < slid::kNumEnergyBin - 1; ig++) {
460  if (fShowerEnergyGev > energyBinBounds[ig] && fShowerEnergyGev < energyBinBounds[ig + 1]) {
461  igev0 = ig;
462  igev1 = ig + 1;
463  }
464  }
465  }
466  if (igev0 > slid::kNumEnergyBin - 1) igev0 = slid::kNumEnergyBin - 1;
467  if (igev1 > slid::kNumEnergyBin - 1) igev1 = slid::kNumEnergyBin - 1;
468  e0 = fShowerEnergyGev - energyBinBounds[igev0];
469  e1 = energyBinBounds[igev1] - fShowerEnergyGev;
470  }
471 
472  //*********************************************************************************************
473  double ParticleIDAlg::CalcTrkHitPath(unsigned int plane)
474  {
475 
476  TVector3 trkDirection = fShower->Dir();
477  TVector3 trkstart = fShower->Start();
478  TVector3 trkstop = fShower->Stop();
479  double trklength = fabs((trkstart - trkstop).Mag());
480  double cellDepth = 2 * fGeom->Plane(plane)->Cell(0)->HalfD();
481 
482  // if track is only one cell in length or short, return the calculated full track length
483  if (fabs(trkDirection.z()) < 1e-6 || fPlaneHits.size() == 1) {
484  return trklength;
485  }
486  // Calculate path length through cell for track with given direction
487  double trkhitpath = cellDepth * (util::pythag(trkDirection.x(),trkDirection.y(),trkDirection.z()) / fabs(trkDirection.z()));
488  return trkhitpath;
489  }
490 
491  //*********************************************************************************************
493  double energySum = 0;
494  std::map<int, int> planecells = fPlaneHits[iPlane];
495  for (std::map<int, int>::iterator itr=planecells.begin(); itr!= planecells.end(); ++itr) {
496  double cellE = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, itr->second, eventID);
497  if ((cellE == -5) || std::isnan(cellE)) cellE = 0.0;
498  energySum += cellE;
499  }
500  double ecor = fNuEEnergyAlg->ECorrMC(energySum);
501  return energySum/ecor;
502  }
503 
504  //NOTE: The input agruement for this function is the actual plane number in the detector
505  //*********************************************************************************************
506  double ParticleIDAlg::CalcPlaneLongitudinalDedx(unsigned int iPlane)
507  {
508  double dedx = 0.0;
509  //
510  // Lazy caching-- only calculate if fDedxLongByPlane[iPlane] not yet set.
511  //
512  std::map< int, double>::iterator planeIter = fDedxLongByPlane.find(iPlane);
513  //std::cout<<"Asked to calculate long dedx for plane "<<iPlane<<std::endl;
514  if (planeIter == fDedxLongByPlane.end()) {
515  double energyInPlane = CalcEnergyInLongitudinalPlane(iPlane);
516  double pathLength = CalcTrkHitPath(iPlane);
517 
518  if (pathLength > 0.0) dedx = energyInPlane / pathLength;
519 
520  fDedxLongByPlane.insert(planeIter, std::pair< int, double>(iPlane, dedx));
521  //std::cout<<"plane idx: "<<iPlane<<" path: "<<pathLength<<" ep: "<<energyInPlane<<" dedx: "<<dedx<<std::endl;
522  }
523  else dedx = planeIter->second;
524  return dedx;
525  }
526 
527  //note that this function is being implemented in the same fasion as RecoJMShower
528  //once we verify it is doing the right thing there may be ways to simplify
529  //*********************************************************************************************
530  double ParticleIDAlg::CalcPlaneTransverseDedx(unsigned int iTransversePlane)
531  {
532  double dedx = -1.0;
533  //
534  // Lazy caching-- only calculate if fDedxTransByPlane[iPlane] not yet set.
535  //
536  std::map< int, double>::iterator transPlaneIter = fDedxTransByPlane.find(iTransversePlane);
537  if (transPlaneIter == fDedxTransByPlane.end()) {
538  //loop over hits in longitudinal planes to calculate transverse plane quantities
539  for(std::map<int, std::map<int, int> >::iterator itr=fPlaneHits.begin(); itr!= fPlaneHits.end(); ++itr){
540  int cplane = itr->first;
541  std::map<int, int> planecells = itr->second;
542  if (planecells.size() == 0) continue;
543  if (fPlaneCentroid[itr->first] == -999.0) continue;
544  double avgPos = fPlaneCentroid[itr->first];
545  int meanCell = (int)avgPos;
546  double wm = avgPos - (double)meanCell;//energy splitting for the centroid cell
547  double wp = 1-wm; //energy splitting for the centroid cell
548  if(meanCell<0)meanCell=0;
549  if(meanCell>(int)fGeom->Plane(cplane)->Ncells()-1)meanCell=(int)fGeom->Plane(cplane)->Ncells()-1;
550  double dir = 1.;
551  if(fabs(fShower->Dir()[2])>0.000001)dir=fabs(fShower->Dir()[2]);
552  int ic = (int)iTransversePlane + 1;
553  double dph = ((double)ic/dir)+wm;//upper limit for cells included in a transverse cell index
554  double dpl = ((double)(ic-1)/dir)+wm;//lower limit for cells included in a transverse cell index
555  //if the mean cell is inside the detector
556  if(meanCell+(int)dph<(int)fGeom->Plane(cplane)->Ncells()){
557  double pfHigh = dph-(int)dph; //fraction of last cell on high end of transverse plane
558  double pfLow = (int)dpl+1-dpl; //fraction of last cell on low end of transverse plane
559  int cellHigh = meanCell + (int)dph;
560  int cellLow = meanCell + (int)dpl;
561  std::map<int, int> cellsInTransPlane;
562  cellsInTransPlane.insert(planecells.lower_bound(cellLow),
563  planecells.upper_bound(cellHigh));
564  for (std::map<int, int>::iterator itrC = cellsInTransPlane.begin(); itrC != cellsInTransPlane.end(); ++itrC){
565  double cellE = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, itrC->second, eventID);
566  if ((cellE == -5) || std::isnan(cellE)) cellE = 0.0;
567  if (itrC->first == cellLow) fTransEnergy[iTransversePlane] += pfLow*cellE;
568  else if (itrC->first == cellHigh) fTransEnergy[iTransversePlane] += pfHigh*cellE;
569  else fTransEnergy[iTransversePlane] += cellE;
570  }
571 
572  fTransPathLength[iTransversePlane] += CalcTrkHitPath(cplane);
573  }
574  double dmh = ((double)ic/dir)+wp;//upper limit for cells included in a transeverse cell index
575  double dml = ((double)(ic-1)/dir)+wp;//lower limit for cells included in a transeverse cell index
576 
577  if(meanCell-(int)dmh>=0){
578  double pfHigh = dmh-(int)dmh; //fraction of last cell on high end of transverse plane
579  double pfLow = (int)dml+1-dml; //fraction of last cell on low end of transverse plane
580  int cellHigh = meanCell - (int)dmh;
581  int cellLow = meanCell - (int)dml;
582  std::map<int, int> cellsInTransPlane;
583  cellsInTransPlane.insert(planecells.lower_bound(cellHigh),
584  planecells.upper_bound(cellLow));
585  for (std::map<int, int>::iterator itrC = cellsInTransPlane.begin(); itrC != cellsInTransPlane.end(); ++itrC){
586  double cellE = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, itrC->second, eventID);
587  if ((cellE == -5) || std::isnan(cellE)) cellE = 0.0;
588  if (itrC->first == cellLow) fTransEnergy[iTransversePlane] += pfLow*cellE;
589  else if (itrC->first == cellHigh) fTransEnergy[iTransversePlane] += pfHigh*cellE;
590  else fTransEnergy[iTransversePlane] += cellE;
591  }
592  fTransPathLength[iTransversePlane] += CalcTrkHitPath(cplane);
593  }
594  }
595  fTransPathLength[iTransversePlane] /= 2.0;
596  if(fTransPathLength[iTransversePlane]>0.00001)dedx = fTransEnergy[iTransversePlane]/fTransPathLength[iTransversePlane];
597  fDedxTransByPlane[iTransversePlane] = dedx;
598  }
599  else dedx = transPlaneIter->second;
600 
601 
602  return dedx;
603  }
604 
605  // Transverse cell Dedx for each plane
606  //*********************************************************************************************
607  double ParticleIDAlg::CalcCellPlaneTransverseDedx(int iTransversePlane, unsigned int iPlane)
608  {
609  double dedx = 0;
610  double transcellE = 0;
611  //loop over hits in longitudinal planes to calculate transverse plane quantities
612  unsigned int gplane = 9999;
613  std::map<unsigned int, unsigned int>::iterator itr0 = fMapIdxToGlobalPlane.find(iPlane);
614  if (itr0 != fMapIdxToGlobalPlane.end()){
615 // std::cout<<iPlane<<" is global "<<"itr0->second "<<itr0->second<<std::endl;
616  gplane = itr0->second;
617  }else{
618 // std::cout<<iPlane<<" "<<"out fMapIdxToGlobalPlane.end() total "<<fMapIdxToGlobalPlane.size()<<std::endl;
619  return dedx;
620  }
621  int pcount=-1;
622  for(std::map<int, std::map<int, int> >::iterator itr=fPlaneHits.begin(); itr!= fPlaneHits.end(); ++itr){
623  pcount++;
624  if((unsigned)itr->first!=gplane)continue;
625  int cplane = itr->first;
626  std::map<int, int> planecells = itr->second;
627  if (planecells.size() == 0) break;
628  if (fPlaneCentroid[itr->first] == -999.0) break;
629 // double avgPos0 = fPlaneCentroid[itr->first];
630  unsigned int planeidx = 0;
631  if (fShower->Dir()[2] >= 0.0){ //forward going shower
632  planeidx = (unsigned int)pcount;
633  }else{
634  planeidx = (unsigned int)(fPlaneHits.size()-pcount-1);
635  }
636  double avgPos = PlaneHitCell(planeidx);
637 // std::cout<<"avg Pos0, hit Pos "<<avgPos0<<" "<<avgPos<<std::endl;
638  int meanCell = (int)avgPos;
639  double wm = avgPos - (double)meanCell;//energy splitting for the centroid cell
640  double wp = 1-wm; //energy splitting for the centroid cell
641  if(meanCell<0)meanCell=0;
642  if(meanCell>(int)fGeom->Plane(cplane)->Ncells()-1)meanCell=(int)fGeom->Plane(cplane)->Ncells()-1;
643  double dir = 1.;
644  if(fabs(fShower->Dir()[2])>0.000001)dir=fabs(fShower->Dir()[2]);
645  int ic = abs((int)iTransversePlane) + 1;
646  double dph = ((double)ic/dir)+wm;//upper limit for cells included in a transverse cell index
647  double dpl = ((double)(ic-1)/dir)+wm;//lower limit for cells included in a transverse cell index
648  //if the mean cell is inside the detector
649 
650  if(iTransversePlane>=0){
651  if(meanCell+(int)dph<(int)fGeom->Plane(cplane)->Ncells()){
652  double pfHigh = dph-(int)dph; //fraction of last cell on high end of transverse plane
653  double pfLow = (int)dpl+1-dpl; //fraction of last cell on low end of transverse plane
654  int cellHigh = meanCell + (int)dph;
655  int cellLow = meanCell + (int)dpl;
656  std::map<int, int> cellsInTransPlane;
657  cellsInTransPlane.insert(planecells.lower_bound(cellLow),
658  planecells.upper_bound(cellHigh));
659  for (std::map<int, int>::iterator itrC = cellsInTransPlane.begin(); itrC != cellsInTransPlane.end(); ++itrC){
660  double cellE = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, itrC->second, eventID, false);
661  if ((cellE == -5) || std::isnan(cellE)) cellE = 0.0;
662  if (itrC->first == cellLow) transcellE += pfLow*cellE;
663  else if (itrC->first == cellHigh) transcellE += pfHigh*cellE;
664  else transcellE += cellE;
665  }
666  }
667  } else {
668 
669  double dmh = ((double)ic/dir)+wp;//upper limit for cells included in a transeverse cell index
670  double dml = ((double)(ic-1)/dir)+wp;//lower limit for cells included in a transeverse cell index
671  if(meanCell-(int)dmh>=0){
672  double pfHigh = dmh-(int)dmh; //fraction of last cell on high end of transverse plane
673  double pfLow = (int)dml+1-dml; //fraction of last cell on low end of transverse plane
674  int cellHigh = meanCell - (int)dmh;
675  int cellLow = meanCell - (int)dml;
676  std::map<int, int> cellsInTransPlane;
677  cellsInTransPlane.insert(planecells.lower_bound(cellHigh),
678  planecells.upper_bound(cellLow));
679  for (std::map<int, int>::iterator itrC = cellsInTransPlane.begin(); itrC != cellsInTransPlane.end(); ++itrC){
680  double cellE = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, itrC->second, eventID, false);
681  if ((cellE == -5) || std::isnan(cellE)) cellE = 0.0;
682  if (itrC->first == cellLow) transcellE += pfLow*cellE;
683  else if (itrC->first == cellHigh) transcellE += pfHigh*cellE;
684  else transcellE += cellE;
685  }
686  }
687  }
688  break;
689  }
690  double path = CalcTrkHitPath(gplane);
691  if(path>0.00001)dedx = transcellE/path;
692  return dedx;
693  }
694 
695  //*********************************************************************************************
696  double ParticleIDAlg::CalcPlaneLongDedxProb(DedxParticleType particleType, int iGev, unsigned int iPlane, unsigned int plane)
697  {
698  double prob = 0;
699 
700  if (iReg < 0 || iReg > kNumXYRegion - 1) return 1;
701  if (iGev >= kNumEnergyBin || iGev < 0) return 1;
702  if (iPlane >= kNumLongitudinalPlane) return 1;
703  double dedx = CalcPlaneLongitudinalDedx(plane);
704  double normHist = double(fDedxDistributions[particleType].GetLongitudinalPlaneHitDedxHist(iReg, iGev, iPlane)->GetEntries());
705  double dedxHistValue = double(fDedxDistributions[particleType].GetLongitudinalPlaneHitDedxHist(iReg, iGev, iPlane)->Interpolate(dedx));
706  double nHistBins = fDedxDistributions[particleType].GetLongitudinalPlaneHitDedxHist(iReg, iGev, iPlane)->GetNbinsX();
707  if (dedxHistValue < 0.0001) dedxHistValue = 0.0001;
708  if (normHist >= 1) prob = dedxHistValue * nHistBins / normHist;
709  return prob;
710  }
711 
712  //*********************************************************************************************
713  double ParticleIDAlg::CalcPlaneTransDedxProb(DedxParticleType particleType, int iGev, unsigned int iTransversePlane)
714  {
715  double prob = 0;
716  if (iReg < 0 || iReg > kNumXYRegion - 1)return 1;
717  if (iGev >= kNumEnergyBin || iGev < 0)return 1;
718  if (iTransversePlane >= kNumTransversePlane)return 1;
719  double dedx = CalcPlaneTransverseDedx(iTransversePlane);
720  double normHist = double(fDedxDistributions[particleType].GetTransversePlaneCellDedxHist(iReg, iGev, iTransversePlane)->GetEntries());
721  double dedxHistValue = double(fDedxDistributions[particleType].GetTransversePlaneCellDedxHist(iReg, iGev, iTransversePlane)->Interpolate(dedx));
722  double nHistBins = fDedxDistributions[particleType].GetTransversePlaneCellDedxHist(iReg, iGev, iTransversePlane)->GetNbinsX();
723  if (dedxHistValue < 0.0001)dedxHistValue = 0.0001;
724  if (normHist >= 1) prob = dedxHistValue * nHistBins / normHist;
725  return prob;
726  }
727 
728  //*********************************************************************************************
729  double ParticleIDAlg::CalcInterPlaneDedxProb(const DedxParticleType type, unsigned int iLongPlane, unsigned int plane)
730  {
731 
732 
733  double probPXPY = 0;
734  double cos = fShower->Dir()[2];
735  unsigned int plane0 = int( int(iLongPlane) / fabs(cos));
736  unsigned int plane1 = int( int(iLongPlane) / fabs(cos)) + 1;
737  double r0 = double(iLongPlane / fabs(cos)) - double(plane0);
738  double r1 = 1 - r0;
739  double probe0r0 = 1;
740  double probe0r1 = 1;
741  double probe1r0 = 1;
742  double probe1r1 = 1;
743  probe0r0 = CalcPlaneLongDedxProb(type, igev0, plane0, plane);
744  probe0r1 = CalcPlaneLongDedxProb(type, igev0, plane1, plane);
745  probe1r0 = CalcPlaneLongDedxProb(type, igev1, plane0, plane);
746  probe1r1 = CalcPlaneLongDedxProb(type, igev1, plane1, plane);
747 
748  double probe0 = (r1 * probe0r0 + r0 * probe0r1);
749  double probe1 = (r1 * probe1r0 + r0 * probe1r1);
750  if (igev0 != igev1)probPXPY = (e1 * probe0 + e0 * probe1) / (e0 + e1);
751  if (igev0 == igev1)probPXPY = probe0;
752  return probPXPY;
753  }
754 
755  //*********************************************************************************************
756  double ParticleIDAlg::CalcInterCellTransDedxProb(const DedxParticleType type, unsigned int iTransPlane)
757  {
758  double probPXPY = 0;
759  double probe0 = CalcPlaneTransDedxProb(type, igev0, iTransPlane);
760  double probe1 = CalcPlaneTransDedxProb(type, igev1, iTransPlane);
761  if (igev0 != igev1)probPXPY = (e1 * probe0 + e0 * probe1) / (e0 + e1);
762  if (igev0 == igev1)probPXPY = probe0;
763  return probPXPY;
764  }
765 
766  //*********************************************************************************************
767  double ParticleIDAlg::GetGapVertexToShowerStart(const rb::Shower* vShower, TVector3 evtVtx, const art::Event& evt)
768  {
769 
770  //initialize event vertex to start of track and first plane in track
771  //TVector3 evtVtx(vShower->Start()[0], vShower->Start()[1], vShower->Start()[2]);
772  double vtxDoca = PointDoca(evtVtx, vShower->Start(), vShower->Stop());
773  double vtxToStart = (evtVtx - vShower->Start()).Mag();
774  double vtxIntToStart = sqrt(vtxToStart * vtxToStart - vtxDoca * vtxDoca);
775  //(Jianming's comments) record the gap from vertex to the start of the track
776  //if the vertex is inside a track record distance of closest approach
777  //TODO: For fuzzyK there can be an internal gap, where the start of the track is misleading
778  //Ruth has looked into this, will try to incorporate internal gaps into this algorithm
779  double shGap = vtxIntToStart;
780  if (vtxIntToStart >= 0 && vtxIntToStart < vtxDoca)shGap = vtxDoca;
781  if (vtxIntToStart < 0)shGap = vtxDoca;
782  return shGap;
783  }
784 
785  //*********************************************************************************************
786  double ParticleIDAlg::PointDoca(TVector3 vtx, TVector3 start, TVector3 stop)
787  {
788  double d = 9999;
789  double denorm = (start - stop).Mag();
790  if (denorm < 0.00001) {
791  d = (vtx - start).Mag();
792  } else {
793  d = ((vtx - start).Cross(vtx - stop)).Mag() / denorm;
794  }
795  return d;
796  }
797 
798  //*********************************************************************************************
799  double ParticleIDAlg::Pi0Mass(const std::vector<const rb::Shower*> showerCollection,
800  unsigned int iShower, const art::Event& evt)
801  {
802  // Load shower if this is a new one.
803  fClosestPi0Mass = -999.;
804  SetShower(showerCollection.at(iShower), showerCollection, evt);
805 
806  double thisShowerERaw = fNuEEnergyAlg->ShowerEnergy(showerCollection[iShower], showerCollection, evt);
807  if ((thisShowerERaw < 1e-10) || std::isnan(thisShowerERaw)) thisShowerERaw = 0.0;
808  TLorentzVector thisShower4Vec(showerCollection.at(iShower)->Dir()[0] * thisShowerERaw,
809  showerCollection.at(iShower)->Dir()[1] * thisShowerERaw,
810  showerCollection.at(iShower)->Dir()[2] * thisShowerERaw, thisShowerERaw);
811 
812  for (unsigned int j = 0; j < showerCollection.size(); j++) { // pi0 mass calculation
813  if (j == iShower)continue;
814  // NOT NECESSARY? if (showerCollection.at(j).ISlice() != showerCollection.at(iShower).ISlice())continue;
815  double showereraw = fNuEEnergyAlg->ShowerEnergy(showerCollection[j], showerCollection, evt);
816  if ((showereraw < 1e-10) || std::isnan(showereraw)) showereraw = 0.0;
817  TLorentzVector otherShower4Vec(showerCollection.at(j)->Dir()[0] * showereraw,
818  showerCollection.at(j)->Dir()[1] * showereraw,
819  showerCollection.at(j)->Dir()[2] * showereraw, showereraw);
820  TLorentzVector p2g = thisShower4Vec + otherShower4Vec;
821  double mgg = p2g.M();
822  if (fabs(fClosestPi0Mass - 0.15) > fabs(mgg - 0.15)) {
823  fClosestPi0Mass = mgg;
824  fPi0Index = j;
825  }
826  }
827  if (fClosestPi0Mass < 0)fClosestPi0Mass = 0.;
828  return fClosestPi0Mass;
829  }
830 
831 
832  ///\brief Index in showerCollection of the shower for which the reconstructed
833  /// invariant mass was closest to pi0 mass for the shower at index iShower.
834  //*********************************************************************************************
835  int ParticleIDAlg::Pi0SecondPhoton(const std::vector< const rb::Shower* > showerCollection,
836  unsigned int iShower, const art::Event& evt)
837  {
838  if( fShower == NULL)
839  Pi0Mass( showerCollection, iShower, evt);
840 
841  return fPi0Index;
842 
843  }// end of Pi0SecondPhoton
844 
845  //*********************************************************************************************
846  double ParticleIDAlg::Radius(const std::vector<const rb::Shower*> showerCollection,
847  unsigned int iShower, const art::Event& evt)
848  {
849  SetShower(showerCollection.at(iShower), showerCollection, evt);
850  double rad = 0;
851  double gev = 0;
852  double dist, di;
853  for (unsigned int nh=0;nh<fShower->NCell();nh++) {
854  art::Ptr<rb::CellHit> cHit = fShower->Cell(nh);
855  double cellEnergy = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, nh, eventID);;
856  TVector3 start = fShower->Start();
857  TVector3 dir = fShower->Dir();
858  TVector3 xyz;
859  fGeom->Plane(cHit->Plane())->Cell(cHit->Cell())->GetCenter(xyz);
860  const TVector3 dxyz = (cHit->View() == geo::kX) ? TVector3(0, 1, 0): TVector3(1, 0, 0);
861  dist = sqrt(geo::ClosestApproach(start, start+dir, xyz, xyz+dxyz, &di));
862  rad += dist*cellEnergy;
863  gev += cellEnergy;
864  }
865  if(gev>0.001)rad = rad/gev;
866  return rad;
867  }
868 
869  //*********************************************************************************************
870  double ParticleIDAlg::PlaneLongDedx(unsigned int pIdx)
871  {
872  std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.find(pIdx);
873  double dedx = -1.0;
874  if (itr != fMapIdxToGlobalPlane.end()) dedx = CalcPlaneLongitudinalDedx(itr->second);
875  else std::cout<<"Requested longitudinal plane index not found, returning dedx=-1"<<std::endl;
876 
877  return dedx;
878  }
879 
880  //*********************************************************************************************
881  double ParticleIDAlg::PlaneTransDedx(unsigned int tpIdx)
882  {
883  double dedx = -1.0;
884  if (tpIdx < slid::kNumTransversePlane) dedx = CalcPlaneTransverseDedx(tpIdx);
885  else std::cout<<"Requested transverse plane index exceeds valid range, maximum plane is: "<<slid::kNumTransversePlane<<", returning dedx=-1"<<std::endl;
886 
887  return dedx;
888  }
889 
890  //*********************************************************************************************
891  double ParticleIDAlg::CellPlaneDedx(int tpIdx, unsigned int pIdx)
892  {
893  double dedx = -1.0;
894  if (tpIdx < slid::kNumTransversePlane) dedx = CalcCellPlaneTransverseDedx(tpIdx, pIdx);
895  else std::cout<<"Requested transverse cell plane index exceeds valid range, maximum plane is: "<<slid::kNumTransversePlane<<", returning dedx=-1"<<std::endl;
896 
897  return dedx;
898  }
899 
900  //*********************************************************************************************
901  double ParticleIDAlg::PlaneRadius(unsigned int pIdx)
902  {
903  std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.find(pIdx);
904  double rad = -1.0;
905  double gev = 0;
906  double dist, di;
907  if (itr != fMapIdxToGlobalPlane.end()){
908  std::map<int, int> planecells = fPlaneHits[itr->second];
909  rad = 0.0;
910  for (std::map<int, int>::iterator it=planecells.begin(); it!= planecells.end(); ++it) {
911  art::Ptr<rb::CellHit> cHit = fShower->Cell(it->second);
912  double cellEnergy = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, it->second, eventID);
913  TVector3 start = fShower->Start();
914  TVector3 dir = fShower->Dir();
915  TVector3 xyz;
916  fGeom->Plane(cHit->Plane())->Cell(cHit->Cell())->GetCenter(xyz);
917  const TVector3 dxyz = (cHit->View() == geo::kX) ? TVector3(0, 1, 0): TVector3(1, 0, 0);
918  dist = sqrt(geo::ClosestApproach(start, start+dir, xyz, xyz+dxyz, &di));
919  rad += dist*cellEnergy;
920  gev += cellEnergy;
921  }
922  }
923  else std::cout<<"Requested longitudinal plane index not found, returning dedx=-1"<<std::endl;
924 
925  if(gev>0.001)rad = rad/gev;
926  return rad;
927  }
928 
929  //*********************************************************************************************
930  unsigned int ParticleIDAlg::PlaneToGlobalPlane(unsigned int pIdx)
931  {
932  std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.find(pIdx);
933  if(itr != fMapIdxToGlobalPlane.end()) return itr->second;
934  return 99999;
935  }
936 
937  //*********************************************************************************************
938  TVector3 ParticleIDAlg::PlaneHitXYZ(unsigned int pIdx)
939  {
940  TVector3 trkplanepos;
941  std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.find(pIdx);
942  if (itr != fMapIdxToGlobalPlane.end()){
943  Double_t cellXYZ[3]={0,0,0};
945  geom->Plane(itr->second)->Cell(0)->GetCenter(cellXYZ);
946  TVector3 trkdir=fShower->Dir();
947  TVector3 trkstart=fShower->Start();
948  TVector3 trkstop=fShower->Stop();
949  double trkplanex=99999,trkplaney=99999, trkplanez=0;
950  if(fabs(trkdir.z())>0.0001){
951  trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
952  trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
953  }
954  else{
955  trkplanex = (trkstart.x()+trkstop.x())/2.0;
956  trkplaney = (trkstart.y()+trkstop.y())/2.0;
957  }
958  trkplanez = cellXYZ[2];
959  trkplanepos.SetX(trkplanex);
960  trkplanepos.SetY(trkplaney);
961  trkplanepos.SetZ(trkplanez);
962  }
963 
964  return trkplanepos;
965  }
966 
967  //*********************************************************************************************
968  double ParticleIDAlg::PlaneHitCell(unsigned int pIdx)
969  {
970  double trkcellId = 9999.0;
971  std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.find(pIdx);
972  if (itr != fMapIdxToGlobalPlane.end()){
973  unsigned int plane = itr->second;
974  Double_t cellXYZ[3]={0,0,0};
975  Double_t cellXYZ1[3];
976  Double_t cellXYZmax[3];
978  geom->Plane(plane)->Cell(0)->GetCenter(cellXYZ);
979  geom->Plane(plane)->Cell(geom->Plane(plane)->Ncells()-1)->GetCenter(cellXYZ1);
980  geom->Plane(plane)->Cell(geom->Plane(plane)->Ncells()-1)->GetCenter(cellXYZmax);
981  double cellW = 2*geom->Plane(plane)->Cell(0)->HalfW();
982  double cellwx = (cellXYZ1[0] - cellXYZ[0] + cellW)/geom->Plane(plane)->Ncells();
983  double cellwy = (cellXYZ1[1] - cellXYZ[1] + cellW)/geom->Plane(plane)->Ncells();
984  TVector3 trkdir=fShower->Dir();
985  TVector3 trkstart=fShower->Start();
986  TVector3 trkstop=fShower->Stop();
987  double trkplanex=99999,trkplaney=99999;//, trkplanez=0;
988  // double trkplanex0=99999,trkplaney0=99999;
989  // double trkplanex1=99999,trkplaney1=99999;
990  if(fabs(trkdir.z())>0.0001){
991  trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
992  trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
993  // trkplanex1 = cellXYZmax[0];
994  // trkplaney1 = cellXYZmax[1];
995  }else{
996  trkplanex = (trkstart.x()+ trkstop.x())/2.0;
997  trkplaney = (trkstart.y()+ trkstop.y())/2.0;
998  }
999  // trkplanez = cellXYZ[2];
1000  geo::View_t view = geom->Plane(plane)->View();
1001  /*
1002  if(view == geo::kX){
1003  std::cout<<"x view plane"<<std::endl;
1004  std::cout<<"x1-x0,cellW "<<" "<<cellXYZ1[0] - cellXYZ[0]<<" "<<cellW<<std::endl;
1005  }else if(view == geo::kY){
1006  std::cout<<"y view plane"<<std::endl;
1007  std::cout<<"y1-y0,cellW "<<" "<<cellXYZ1[1] - cellXYZ[1]<<" "<<cellW<<std::endl;
1008  }
1009  */
1010  // int trkcellId0 = -9999;
1011  // int trkcellId1 = -9999;
1012  // int trkcellIdmax = geom->Plane(plane)->Ncells()-1;
1013  if(view == geo::kX)trkcellId = (trkplanex+cellwx/2-cellXYZ[0])/cellwx;
1014  if(view == geo::kY)trkcellId = (trkplaney+cellwy/2-cellXYZ[1])/cellwy;
1015  // if(view == geo::kX)trkcellId0 = int((trkplanex0+cellwx/2-cellXYZ[0])/cellwx);
1016  // if(view == geo::kY)trkcellId0 = int((trkplaney0+cellwy/2-cellXYZ[1])/cellwy);
1017  // if(view == geo::kX)trkcellId1 = int((trkplanex1+cellwx/2-cellXYZ[0])/cellwx);
1018  // if(view == geo::kY)trkcellId1 = int((trkplaney1+cellwy/2-cellXYZ[1])/cellwy);
1019  }else{std::cout<<pIdx<<" end !!! "<<std::endl;}
1020  return trkcellId;
1021  }
1022 
1023  //*********************************************************************************************
1024  double ParticleIDAlg::PlaneCentroid(unsigned int pIdx)
1025  {
1026  double centroid = 9999.0;
1027  std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.find(pIdx);
1028  if (itr != fMapIdxToGlobalPlane.end()){
1029  std::map<int, double>::iterator it = fPlaneCentroid.find(itr->second);
1030  if (it != fPlaneCentroid.end()) centroid = it->second;
1031  }
1032 
1033  return centroid;
1034  }
1035 
1036  //*********************************************************************************************
1038  {
1040  double I11 = 0;
1041  double I12 = 0;
1042  double I21 = 0;
1043  double I22 = 0;
1044 
1045  // Using 3-D nodes to calculate shower balance
1046  // TVector3 core(shower.Start()[0],shower.Start()[1],shower.Start()[2]); // original point of shower frame
1047  TVector3 iz = fShower->Dir().Unit(); // z in shower frame
1048  TVector3 ix0(1.0,0.0,0);
1049  TVector3 iy0(0.0,1.0,0);
1050  TVector3 ix = iy0.Cross(iz).Unit(); // x axis in shower frame
1051  TVector3 iy = iz.Cross(ix).Unit(); // y axis in shower frame
1052 
1053  std::map<unsigned int, unsigned int>::iterator itr2 = fMapIdxToGlobalPlane.begin();
1054  for(std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.begin(); itr != fMapIdxToGlobalPlane.end(); ++itr){
1055  itr2++;
1056  std::map<int, int> planecells1 = fPlaneHits[itr->second];
1057  if (itr2 == fMapIdxToGlobalPlane.end()) break;
1058  std::map<int, int> planecells2 = fPlaneHits[itr2->second];
1059  if(planecells1.size()==0)continue;
1060  unsigned int plane1 = itr->second;
1061  if(planecells2.size()==0)continue;
1062  unsigned int plane2 = itr2->second;
1063  geo::View_t view1 = geom->Plane(plane1)->View();
1064  geo::View_t view2 = geom->Plane(plane2)->View();
1065  if(view1==view2)continue;
1066  TVector3 core = PlaneHitXYZ(itr->first);
1067  TVector3 pcore = PlaneHitXYZ(itr2->first);
1068  std::vector<TVector3> planenodecol;
1069  std::vector<double> planenodeenergycol;
1070  planenodecol.clear();
1071  planenodeenergycol.clear();
1072  for (std::map<int, int>::iterator citr = planecells1.begin(); citr != planecells1.end(); ++citr) {
1073  unsigned int cell1 = citr->first;
1074  for(std::map<int, int>::iterator citr2 = planecells2.begin(); citr2 != planecells2.end(); ++citr2){
1075  unsigned int cell2 = citr2->first;
1076  TVector3 node = GetCellNodePos(plane1, cell1, plane2, cell2);// define 3-D node
1077  double nodez = (pcore - fShower->Start()).Mag();
1078  TVector3 nodetopcore = node - pcore;
1079  nodetopcore.SetZ(nodez);
1080  TVector3 nodetocore = node - core;
1081  nodetocore.SetZ(nodez);
1082 
1083  TVector3 xyz = node - core;//transform xy to shower frame
1084 
1085  double p = fShower->Dir()[0];
1086  double q = fShower->Dir()[1];
1087  double r = fShower->Dir()[2];
1088  double x0 = xyz[0];
1089  double y0 = xyz[1];
1090  double z0 = xyz[2];
1091  TVector3 vt(q*z0-r*y0,r*x0-p*z0,p*y0-q*x0);
1092 
1093  double nodeenergy = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, citr->second, eventID) +
1094  fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, citr2->second, eventID);;
1095 
1096  I11 += nodeenergy*(nodetocore.Y()*nodetocore.Y());
1097  I12 += -nodeenergy*(nodetocore.X()*nodetocore.Y());
1098  I21 += -nodeenergy*(nodetocore.X()*nodetocore.Y());
1099  I22 += nodeenergy*(nodetocore.X()*nodetocore.X());
1100 
1101  }
1102  }
1103  }
1104  // calculate eiganvalues and diagonalize inertial tensor
1105  double Ixx = 0.5*(I11+I22+sqrt(pow(I11+I22,2)-4*(I11*I22-I12*I21)));
1106  double Iyy = 0.5*(I11+I22-sqrt(pow(I11+I22,2)-4*(I11*I22-I12*I21)));
1107  double asym = 1;
1108  if(fabs(Ixx+Iyy)>1e-10)asym=(Ixx-Iyy)/(Ixx+Iyy);
1109 
1110  TVector3 ixyz(Ixx,Iyy,0);
1111  fIxyz = ixyz;
1112  fAsym = asym;
1113  }
1114 
1115  //*********************************************************************************************
1116  TVector3 ParticleIDAlg::GetCellNodePos(unsigned int plane1, unsigned int cell1, unsigned int plane2, unsigned int cell2)
1117  {
1118  double cellXYZ1[3]={0,0,0};
1119  double cellXYZ2[3]={0,0,0};
1121  geom->Plane(plane1)->Cell(cell1)->GetCenter(cellXYZ1);
1122  geom->Plane(plane2)->Cell(cell2)->GetCenter(cellXYZ2);
1123  geo::View_t view1 = geom->Plane(plane1)->View();
1124  geo::View_t view2 = geom->Plane(plane2)->View();
1125  TVector3 nullpos(-9999,-9999,-9999);
1126  if(view1==view2)return nullpos;
1127  double x=0,y=0,z=0;
1128  if(view1 == geo::kX && view2 == geo::kY){
1129  x = cellXYZ1[0];
1130  y = cellXYZ2[1];
1131  }else if(view1 == geo::kY && view2 == geo::kX){
1132  x = cellXYZ2[0];
1133  y = cellXYZ1[1];
1134  }
1135  z = (cellXYZ1[2]+cellXYZ2[2])/2.;
1136  TVector3 nodepos(x,y,z);
1137  return nodepos;
1138  }
1139 
1140  //*********************************************************************************************
1142  const std::vector< const rb::Shower* > showercol,
1143  const art::Event& evt)
1144  {
1145  if(shower == NULL){
1146  return 0;
1147  }
1148  SetShower( shower, showercol, evt);
1149 
1150  int nMipPlanes = 0;
1151  int nPlanes = 20;
1152  if( fPlaneHits.size() < 20 )
1153  nPlanes = fPlaneHits.size();
1154  if(nPlanes == 0)
1155  return 0;
1156 
1157  int count = 0;
1158 
1159  for( auto const & iPlaneHit : fPlaneHits ){
1160  if(count > nPlanes)
1161  break;
1162  ++count;
1163 
1164  double dedx = CalcPlaneLongitudinalDedx( iPlaneHit.first );
1165 
1166  if( dedx > fMipRangeHigh || dedx < fMipRangeLow )
1167  continue;
1168 
1169  nMipPlanes++;
1170  }
1171  return nMipPlanes;
1172  }// end of GetMIPPlanes
1173 
1174 }// End of namespace slid
novadaq::cnv::DetId fDetId
Detector ID.
NuEEnergyAlg * fNuEEnergyAlg
Nue energy algorithms.
static const int kNumTransversePlane
Number of "transverse" planes considered for dE/dx histograms ("planes" transverse to shower axis)...
double HalfD() const
Definition: CellGeo.cxx:205
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
std::map< int, double > fTransEnergy
Map of transverse energy for each transverse plane.
double CalcPlaneTransDedxProb(DedxParticleType particleType, int igev, unsigned int iTransversePlane)
Get transverse dE/dx probability for a given transverse plane, energy bin, and particle hypothesis...
std::map< int, std::map< int, int > > CalcPlaneHits(const rb::Shower *vShower)
Return a map off all hits in a given longitudinal plane indexed by the longitudinalplane for shower a...
set< int >::iterator it
double fClosestPi0Mass
Closest mass to pi0 mass formed by this shower and any other shower in the slice. ...
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
int iReg
Detector sub-region bin. Calculated once and cached for a given shower/event combo.
double CalcPlaneTransverseDedx(unsigned int iTransversePlane)
Calculate the transverse dE/dx in a given tranverse slice. Input is the transverse slice index...
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
double DedxTransLL(const slid::DedxParticleType partHyp, const rb::Shower *vShower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Calculate transverse dE/dx log-likelihood for specific particle hypothesis.
unsigned short Plane() const
Definition: CellHit.h:39
double GetGapVertexToShowerStart(const rb::Shower *vShower, TVector3 evtVtx, const art::Event &evt)
Calculate the gap between the vertex and the start of the shower.
geo::View_t View() const
Definition: CellHit.h:41
TVector3 PlaneHitXYZ(unsigned int pIdx)
Return point of intersection between shower core and plane.
double CalcCellPlaneTransverseDedx(int iTransversePlane, unsigned int iPlane)
Calculate the transverse dE/dx in a given tranverse cell and plane. Input is the transverse slice ind...
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
bool equal(T *first, T *second)
double e0
Energy values of energy bins adjacent to shower. Calculated once and cached for a given shower/event ...
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
void CalcAsymIneria()
Calculate the asymmetry and inertia of the shower.
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
double HalfW() const
Definition: CellGeo.cxx:191
double PlaneRadius(unsigned int pIdx)
return shower radius for a plane index
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void abs(TH1 *hist)
double CalcInterPlaneDedxProb(DedxParticleType type, unsigned int iLongPlane, unsigned int plane)
From Jianming&#39;s code, still working out changes.
static const int kNumXYRegion
Number of XY regions into which detector is divided for dE/dx histograms.
std::map< unsigned int, unsigned int > fMapIdxToGlobalPlane
std::vector< slid::DedxDistribution > fDedxDistributions
dE/dx distributions (histogram files and data) for each detector type and particle category ...
double fShowerEnergyGev
Shower energy.
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
double PlaneLongDedx(unsigned int pIdx)
return longitudinal dedx for a specified plane index
cout<< t1-> GetEntries()
Definition: plottest35.C:29
virtual void SetStart(TVector3 start)
Definition: Prong.cxx:75
void reconfigure(const fhicl::ParameterSet &pset, novadaq::cnv::DetId detId)
Reinitialize parameters read from FCL file.
double fMipRangeLow
Lower bound of MIP Range.
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
TVector3 GetCellNodePos(unsigned int plane1, unsigned int cell1, unsigned int plane2, unsigned int cell2)
Helper function for interia and asymmetry calculation.
double dist
Definition: runWimpSim.h:113
int GetMIPPlanes(const rb::Shower *shower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Number of MIP planes within the first 20 planes of the start of the shower.
unsigned int PlaneToGlobalPlane(unsigned int pIdx)
Return the gloabel plane index for an index in the shower coordinate.
unsigned short Cell() const
Definition: CellHit.h:40
DedxParticleType
An enum used to give allowed particle types a visible name in the code. Note that for electron...
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
std::vector< const rb::Shower * > fShowerCol
double CellPlaneDedx(int tpIdx, unsigned int pIdx)
return transverse dedx for specified transverse plane
int Pi0SecondPhoton(const std::vector< const rb::Shower * > showerCollection, unsigned int iShower, const art::Event &evt)
Index in showerCollection of the shower for which the reconstructed invariant mass was closest to pi0...
double PlaneCentroid(unsigned int pIdx)
Return the energy weighted mean shower positon for the plane in cm.
Prototype Near Detector on the surface at FNAL.
double CalcPlaneLongitudinalDedx(unsigned int plane)
Calculate the longitudinal dE/dx in a given plane. Input is the actuall plane number in the detector...
virtual void SetDir(TVector3 dir)
Definition: Prong.cxx:104
T get(std::string const &key) const
Definition: ParameterSet.h:231
double CellEnergy(const rb::Shower *shower, const std::vector< const rb::Shower * > showercol, const int &index, const art::EventID &evtid, bool useweight=true, bool *isCalibrated=NULL)
Returns the energy deposited in a cell. Used in computation of total shower energy. Only used if fUseStdCellE is false.
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
TVector3 Stop() const
Endpoint of the shower.
Definition: Shower.cxx:55
int evt
std::map< int, double > CalcPlaneCentroid(const rb::Shower *vShower)
Return a map of mean cell posisiton for each plane for shower already read into object by SetShower...
ParticleIDAlg(fhicl::ParameterSet const &pset, novadaq::cnv::DetId detId, NuEEnergyAlg *NuEEnergyAlgIn, bool loadLibs=true)
bool LoadDedxHistogramFiles(const fhicl::ParameterSet &pset, novadaq::cnv::DetId geom)
Load all of the histograms from files into arrays of TH1D* objects. The details of handling the histo...
Collect Geo headers and supply basic geometry functions.
Float_t d
Definition: plot.C:236
TVector3 pos
Shower start position.
void cellPos()
Definition: cellShifts.C:26
def Interpolate(x1, y1, x2, y2, yvalue)
Definition: HandyFuncs.py:218
const double j
Definition: BetheBloch.cxx:29
art::ServiceHandle< geo::Geometry > fGeom
geometry
std::map< int, float > fPartTransLL
Map of the transverse ll by particle type.
static constexpr Double_t rad
Definition: Munits.h:162
void CalculateEnergyBinIndex()
Calculate energy bin indices and energy values from the histograms corresponding energy of shower to ...
double PointDoca(TVector3 vtx, TVector3 start, TVector3 stop)
Get distance of closest approach between shower and vertex.
z
Definition: test.py:28
double fAsym
shower asymmetry
art::EventID eventID
data/MC flag
OStream cout
Definition: OStream.cxx:6
double ShowerEnergy(const rb::Shower *shower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Returns the total energy of a shower. along with corrections due to dead material and threshold effec...
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
const std::string path
Definition: plot_BEN.C:43
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
double PlaneTransDedx(unsigned int tpIdx)
return transverse dedx for specified transverse plane
std::map< int, double > fTransPathLength
Map of path length for each transverse plane.
Calculate dE/dx and log likelihoods for different particle hypotheses. This information will be of us...
std::map< int, double > fDedxTransByPlane
double CalcEnergyInLongitudinalPlane(int iPlane)
Sum the energy in a longitudinal plane.
double ClosestApproach(const double point[], const double intercept[], const double slopes[], double closest[])
Find the distance of closest approach between point and line.
Definition: Geo.cxx:222
TVector3 fIxyz
shower moment of ineria
TDirectory * dir
Definition: macro.C:5
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
A rb::Prong with a length.
Definition: Shower.h:18
static const int kNumEnergyBin
Number of energy bins into which detector is divided for dE/dx histograms.
void geom(int which=0)
Definition: geom.C:163
float Mag() const
TFile * LoadHistFromFile(TString filePath)
T cos(T number)
Definition: d0nt_math.hpp:78
double ECorrMC(double gev)
EM shower energy correction for MC Only used if fCorAbsCell is false.
TRandom3 r(0)
int ncells
Definition: geom.C:124
static const float kEnergyBins[]
bin boundaries of energy
double Radius(const std::vector< const rb::Shower * > showerCollection, unsigned int iShower, const art::Event &evt)
Calculate the shower radius.
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
Build slid::LID objects to store electron ID, if asked for, otherwise, calculate LID info and make av...
Definition: FillPIDs.h:13
const int ID() const
Definition: Cluster.h:75
std::map< int, double > fPlaneCentroid
static const int kNumLongitudinalPlane
Number of longitudinal planes considered for dE/dx histograms.
std::map< int, std::map< int, int > > fPlaneHits
double CalcInterCellTransDedxProb(const DedxParticleType type, unsigned int iTransPlane)
From Jianming&#39;s code, still working out changes.
Float_t e
Definition: plot.C:35
std::map< int, float > fPartLongLL
Map of the longitudinal ll by paricle type.
double PlaneHitCell(unsigned int pIdx)
Return point of intersection cell between shower core and plane.
double DedxInverseLongLL(const slid::DedxParticleType partHyp, const rb::Shower *vShower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis assuming that the partic...
double Pi0Mass(const std::vector< const rb::Shower * > showerCollection, unsigned int iShower, const art::Event &evt)
Calculate the mass of candidate and any other shower in the slice that is closest to the pi0 mass...
int igev0
Shower adjacent energy bins. Calculated once and cached for a given shower/event combo.
EventID id() const
Definition: Event.h:56
std::map< int, double > fDedxLongByPlane
void CalcDetectorXYRegionIndex()
Determine detector XY region.
double CalcPlaneLongDedxProb(DedxParticleType particleType, int igev, unsigned int iPlane, unsigned int plane)
Get longitudinal dE/dx probability for a given longitudinal plane, energy bin, and particle hypothesi...
double DedxLongLL(const slid::DedxParticleType partHyp, const rb::Shower *vShower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis.
double CalcTrkHitPath(unsigned int plane)
Calculate path length of shower through this plane.
double fMipRangeHigh
Upper bound of MIP Range.
void SetShower(const rb::Shower *vShower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Set rb::Prong to be analyzed. This must be set before any calculations can be done.
enum BeamMode string