Public Member Functions | Public Attributes | Private Member Functions | Private Attributes | List of all members
slid::ParticleIDAlg Class Reference

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-10-23/ShowerLID/ParticleIDAlg.h"

Public Member Functions

 ParticleIDAlg (fhicl::ParameterSet const &pset, novadaq::cnv::DetId detId, NuEEnergyAlg *NuEEnergyAlgIn, bool loadLibs=true)
 
 ~ParticleIDAlg ()
 
void reconfigure (const fhicl::ParameterSet &pset, novadaq::cnv::DetId detId)
 Reinitialize parameters read from FCL file. More...
 
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. More...
 
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. More...
 
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 particle direction if opposite to the reconstructed direction. More...
 
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. More...
 
double GetGapVertexToShowerStart (const rb::Shower *vShower, TVector3 evtVtx, const art::Event &evt)
 Calculate the gap between the vertex and the start of the shower. More...
 
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. More...
 
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 mass for the shower at index iShower. More...
 
double Radius (const std::vector< const rb::Shower * > showerCollection, unsigned int iShower, const art::Event &evt)
 Calculate the shower radius. More...
 
double PlaneLongDedx (unsigned int pIdx)
 return longitudinal dedx for a specified plane index More...
 
double PlaneTransDedx (unsigned int tpIdx)
 return transverse dedx for specified transverse plane More...
 
double CellPlaneDedx (int tpIdx, unsigned int pIdx)
 return transverse dedx for specified transverse plane More...
 
double PlaneRadius (unsigned int pIdx)
 return shower radius for a plane index More...
 
TVector3 PlaneHitXYZ (unsigned int pIdx)
 Return point of intersection between shower core and plane. More...
 
double PlaneHitCell (unsigned int pIdx)
 Return point of intersection cell between shower core and plane. More...
 
unsigned int PlaneToGlobalPlane (unsigned int pIdx)
 Return the gloabel plane index for an index in the shower coordinate. More...
 
double PlaneCentroid (unsigned int pIdx)
 Return the energy weighted mean shower positon for the plane in cm. More...
 
double ShowerAsym () const
 Return the asymmetry of the shower. More...
 
TVector3 InertiaXYZ () const
 Return the inertia calculation. More...
 
double PointDoca (TVector3 vtx, TVector3 start, TVector3 stop)
 Get distance of closest approach between shower and vertex. More...
 
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. More...
 

Public Attributes

std::map< int, float > fPartLongLL
 Map of the longitudinal ll by paricle type. More...
 
std::map< int, float > fPartTransLL
 Map of the transverse ll by particle type. More...
 

Private Member Functions

double CalcPlaneLongitudinalDedx (unsigned int plane)
 Calculate the longitudinal dE/dx in a given plane. Input is the actuall plane number in the detector. More...
 
double CalcPlaneTransverseDedx (unsigned int iTransversePlane)
 Calculate the transverse dE/dx in a given tranverse slice. Input is the transverse slice index. More...
 
double CalcCellPlaneTransverseDedx (int iTransversePlane, unsigned int iPlane)
 Calculate the transverse dE/dx in a given tranverse cell and plane. Input is the transverse slice index, including plus and minus to the shower core, and plane number. More...
 
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 already read into object by SetShower. More...
 
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. More...
 
double CalcEnergyInLongitudinalPlane (int iPlane)
 Sum the energy in a longitudinal plane. More...
 
void CalcDetectorXYRegionIndex ()
 Determine detector XY region. More...
 
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 hypothesis. More...
 
double CalcPlaneTransDedxProb (DedxParticleType particleType, int igev, unsigned int iTransversePlane)
 Get transverse dE/dx probability for a given transverse plane, energy bin, and particle hypothesis. More...
 
double CalcInterPlaneDedxProb (DedxParticleType type, unsigned int iLongPlane, unsigned int plane)
 From Jianming's code, still working out changes. More...
 
double CalcInterCellTransDedxProb (const DedxParticleType type, unsigned int iTransPlane)
 From Jianming's code, still working out changes. More...
 
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 histograms are taken care of within object of the helper class slid::DedxDistribution. More...
 
double CalcTrkHitPath (unsigned int plane)
 Calculate path length of shower through this plane. More...
 
void CalculateEnergyBinIndex ()
 Calculate energy bin indices and energy values from the histograms corresponding energy of shower to be used to interpolate between bins in energy. More...
 
void CalcAsymIneria ()
 Calculate the asymmetry and inertia of the shower. More...
 
TVector3 GetCellNodePos (unsigned int plane1, unsigned int cell1, unsigned int plane2, unsigned int cell2)
 Helper function for interia and asymmetry calculation. More...
 

Private Attributes

art::EventID eventID
 data/MC flag More...
 
double fClosestPi0Mass
 Closest mass to pi0 mass formed by this shower and any other shower in the slice. More...
 
std::vector< slid::DedxDistributionfDedxDistributions
 dE/dx distributions (histogram files and data) for each detector type and particle category More...
 
const rb::ShowerfShower
 rb::Shower object for which particle ID is to be calculated. This is cached. More...
 
std::vector< const rb::Shower * > fShowerCol
 
double fShowerEnergyGev
 Shower energy. More...
 
art::ServiceHandle< geo::GeometryfGeom
 geometry More...
 
NuEEnergyAlgfNuEEnergyAlg
 Nue energy algorithms. More...
 
std::map< int, double > fDedxLongByPlane
 
std::map< int, double > fDedxTransByPlane
 
std::map< int, std::map< int, int > > fPlaneHits
 
std::map< unsigned int, unsigned intfMapIdxToGlobalPlane
 
std::map< int, double > fPlaneCentroid
 
std::map< int, double > fTransEnergy
 Map of transverse energy for each transverse plane. More...
 
std::map< int, double > fTransPathLength
 Map of path length for each transverse plane. More...
 
int fPi0Index
 
int iReg
 Detector sub-region bin. Calculated once and cached for a given shower/event combo. More...
 
TVector3 pos
 Shower start position. More...
 
int igev0
 Shower adjacent energy bins. Calculated once and cached for a given shower/event combo. More...
 
int igev1
 
double e0
 Energy values of energy bins adjacent to shower. Calculated once and cached for a given shower/event combo. More...
 
double e1
 
TVector3 fIxyz
 shower moment of ineria More...
 
double fAsym
 shower asymmetry More...
 
double fMipRangeHigh
 Upper bound of MIP Range. More...
 
double fMipRangeLow
 Lower bound of MIP Range. More...
 
novadaq::cnv::DetId fDetId
 Detector ID. More...
 

Detailed Description

Definition at line 50 of file ParticleIDAlg.h.

Constructor & Destructor Documentation

slid::ParticleIDAlg::ParticleIDAlg ( fhicl::ParameterSet const &  pset,
novadaq::cnv::DetId  detId,
NuEEnergyAlg NuEEnergyAlgIn,
bool  loadLibs = true 
)
explicit

Definition at line 30 of file ParticleIDAlg.cxx.

References LoadDedxHistogramFiles().

34  : fClosestPi0Mass(-9999),
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){
45  }
46  }
novadaq::cnv::DetId fDetId
Detector ID.
NuEEnergyAlg * fNuEEnergyAlg
Nue energy algorithms.
double fClosestPi0Mass
Closest mass to pi0 mass formed by this shower and any other shower in the slice. ...
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
std::vector< slid::DedxDistribution > fDedxDistributions
dE/dx distributions (histogram files and data) for each detector type and particle category ...
double fMipRangeLow
Lower bound of MIP Range.
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...
double fAsym
shower asymmetry
double fMipRangeHigh
Upper bound of MIP Range.
slid::ParticleIDAlg::~ParticleIDAlg ( )

Definition at line 49 of file ParticleIDAlg.cxx.

50  {
51  }

Member Function Documentation

void slid::ParticleIDAlg::CalcAsymIneria ( )
private

Calculate the asymmetry and inertia of the shower.

Definition at line 1037 of file ParticleIDAlg.cxx.

References slid::NuEEnergyAlg::CellEnergy(), rb::Prong::Dir(), e, eventID, stan::math::fabs(), fAsym, fIxyz, fMapIdxToGlobalPlane, fNuEEnergyAlg, fPlaneHits, fShower, fShowerCol, geom(), GetCellNodePos(), Mag(), file_size_ana::node, geo::GeometryBase::Plane(), PlaneHitXYZ(), cet::pow(), r(), std::sqrt(), rb::Prong::Start(), and geo::PlaneGeo::View().

Referenced by SetShower().

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  }
NuEEnergyAlg * fNuEEnergyAlg
Nue energy algorithms.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
TVector3 PlaneHitXYZ(unsigned int pIdx)
Return point of intersection between shower core and plane.
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
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
std::map< unsigned int, unsigned int > fMapIdxToGlobalPlane
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
TVector3 GetCellNodePos(unsigned int plane1, unsigned int cell1, unsigned int plane2, unsigned int cell2)
Helper function for interia and asymmetry calculation.
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
std::vector< const rb::Shower * > fShowerCol
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
double fAsym
shower asymmetry
art::EventID eventID
data/MC flag
TVector3 fIxyz
shower moment of ineria
void geom(int which=0)
Definition: geom.C:163
float Mag() const
TRandom3 r(0)
std::map< int, std::map< int, int > > fPlaneHits
Float_t e
Definition: plot.C:35
double slid::ParticleIDAlg::CalcCellPlaneTransverseDedx ( int  iTransversePlane,
unsigned int  iPlane 
)
private

Calculate the transverse dE/dx in a given tranverse cell and plane. Input is the transverse slice index, including plus and minus to the shower core, and plane number.

Definition at line 607 of file ParticleIDAlg.cxx.

References abs(), CalcTrkHitPath(), slid::NuEEnergyAlg::CellEnergy(), dir, rb::Prong::Dir(), eventID, stan::math::fabs(), fGeom, fMapIdxToGlobalPlane, fNuEEnergyAlg, fPlaneCentroid, fPlaneHits, fShower, fShowerCol, makeTrainCVSamples::int, std::isnan(), geo::PlaneGeo::Ncells(), path, geo::GeometryBase::Plane(), and PlaneHitCell().

Referenced by CellPlaneDedx().

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  }
NuEEnergyAlg * fNuEEnergyAlg
Nue energy algorithms.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
void abs(TH1 *hist)
std::map< unsigned int, unsigned int > fMapIdxToGlobalPlane
const PlaneGeo * Plane(unsigned int i) const
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
std::vector< const rb::Shower * > fShowerCol
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
art::ServiceHandle< geo::Geometry > fGeom
geometry
art::EventID eventID
data/MC flag
const std::string path
Definition: plot_BEN.C:43
TDirectory * dir
Definition: macro.C:5
std::map< int, double > fPlaneCentroid
std::map< int, std::map< int, int > > fPlaneHits
double PlaneHitCell(unsigned int pIdx)
Return point of intersection cell between shower core and plane.
double CalcTrkHitPath(unsigned int plane)
Calculate path length of shower through this plane.
void slid::ParticleIDAlg::CalcDetectorXYRegionIndex ( )
private

Determine detector XY region.

Definition at line 416 of file ParticleIDAlg.cxx.

References fDetId, fShower, iReg, cmf::kFARDET, pos, rb::Prong::Start(), and rb::Shower::Stop().

Referenced by SetShower().

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  }
novadaq::cnv::DetId fDetId
Detector ID.
int iReg
Detector sub-region bin. Calculated once and cached for a given shower/event combo.
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
virtual TVector3 Start() const
Definition: Prong.h:73
TVector3 Stop() const
Endpoint of the shower.
Definition: Shower.cxx:55
TVector3 pos
Shower start position.
double slid::ParticleIDAlg::CalcEnergyInLongitudinalPlane ( int  iPlane)
private

Sum the energy in a longitudinal plane.

Definition at line 492 of file ParticleIDAlg.cxx.

References slid::NuEEnergyAlg::CellEnergy(), slid::NuEEnergyAlg::ECorrMC(), eventID, fNuEEnergyAlg, fPlaneHits, fShower, fShowerCol, and std::isnan().

Referenced by CalcPlaneLongitudinalDedx().

492  {
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  }
NuEEnergyAlg * fNuEEnergyAlg
Nue energy algorithms.
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
std::vector< const rb::Shower * > fShowerCol
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.
art::EventID eventID
data/MC flag
double ECorrMC(double gev)
EM shower energy correction for MC Only used if fCorAbsCell is false.
std::map< int, std::map< int, int > > fPlaneHits
double slid::ParticleIDAlg::CalcInterCellTransDedxProb ( const DedxParticleType  type,
unsigned int  iTransPlane 
)
private

From Jianming's code, still working out changes.

Definition at line 756 of file ParticleIDAlg.cxx.

References CalcPlaneTransDedxProb(), e0, e1, igev0, and igev1.

Referenced by DedxTransLL().

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  }
double CalcPlaneTransDedxProb(DedxParticleType particleType, int igev, unsigned int iTransversePlane)
Get transverse dE/dx probability for a given transverse plane, energy bin, and particle hypothesis...
double e0
Energy values of energy bins adjacent to shower. Calculated once and cached for a given shower/event ...
int igev0
Shower adjacent energy bins. Calculated once and cached for a given shower/event combo.
double slid::ParticleIDAlg::CalcInterPlaneDedxProb ( DedxParticleType  type,
unsigned int  iLongPlane,
unsigned int  plane 
)
private

From Jianming's code, still working out changes.

Definition at line 729 of file ParticleIDAlg.cxx.

References CalcPlaneLongDedxProb(), std::cos(), rb::Prong::Dir(), e0, e1, stan::math::fabs(), fShower, igev0, igev1, and makeTrainCVSamples::int.

Referenced by DedxLongLL().

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  }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double e0
Energy values of energy bins adjacent to shower. Calculated once and cached for a given shower/event ...
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
T cos(T number)
Definition: d0nt_math.hpp:78
int igev0
Shower adjacent energy bins. Calculated once and cached for a given shower/event combo.
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...
std::map< int, double > slid::ParticleIDAlg::CalcPlaneCentroid ( const rb::Shower vShower)
private

Return a map of mean cell posisiton for each plane for shower already read into object by SetShower.

Definition at line 278 of file ParticleIDAlg.cxx.

References slid::NuEEnergyAlg::CellEnergy(), cellPos(), eventID, fNuEEnergyAlg, fPlaneHits, fShowerCol, and std::isnan().

Referenced by SetShower().

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  }
NuEEnergyAlg * fNuEEnergyAlg
Nue energy algorithms.
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
std::vector< const rb::Shower * > fShowerCol
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.
void cellPos()
Definition: cellShifts.C:26
art::EventID eventID
data/MC flag
std::map< int, std::map< int, int > > fPlaneHits
std::map< int, std::map< int, int > > slid::ParticleIDAlg::CalcPlaneHits ( const rb::Shower vShower)
private

Return a map off all hits in a given longitudinal plane indexed by the longitudinalplane for shower already read into object by SetShower.

Definition at line 241 of file ParticleIDAlg.cxx.

References rb::CellHit::Cell(), rb::Cluster::Cell(), rb::Prong::Dir(), fMapIdxToGlobalPlane, MECModelEnuComparisons::i, rb::Cluster::MaxPlane(), rb::Cluster::MinPlane(), rb::Cluster::NCell(), ncells, and rb::CellHit::Plane().

Referenced by SetShower().

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
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
unsigned short Plane() const
Definition: CellHit.h:39
std::map< unsigned int, unsigned int > fMapIdxToGlobalPlane
unsigned short Cell() const
Definition: CellHit.h:40
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
int ncells
Definition: geom.C:124
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
double slid::ParticleIDAlg::CalcPlaneLongDedxProb ( DedxParticleType  particleType,
int  igev,
unsigned int  iPlane,
unsigned int  plane 
)
private

Get longitudinal dE/dx probability for a given longitudinal plane, energy bin, and particle hypothesis.

Definition at line 696 of file ParticleIDAlg.cxx.

References CalcPlaneLongitudinalDedx(), fDedxDistributions, GetEntries(), HandyFuncs::Interpolate(), iReg, slid::kNumEnergyBin, slid::kNumLongitudinalPlane, and slid::kNumXYRegion.

Referenced by CalcInterPlaneDedxProb().

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  }
int iReg
Detector sub-region bin. Calculated once and cached for a given shower/event combo.
static const int kNumXYRegion
Number of XY regions into which detector is divided for dE/dx histograms.
std::vector< slid::DedxDistribution > fDedxDistributions
dE/dx distributions (histogram files and data) for each detector type and particle category ...
cout<< t1-> GetEntries()
Definition: plottest35.C:29
double CalcPlaneLongitudinalDedx(unsigned int plane)
Calculate the longitudinal dE/dx in a given plane. Input is the actuall plane number in the detector...
def Interpolate(x1, y1, x2, y2, yvalue)
Definition: HandyFuncs.py:218
static const int kNumEnergyBin
Number of energy bins into which detector is divided for dE/dx histograms.
static const int kNumLongitudinalPlane
Number of longitudinal planes considered for dE/dx histograms.
double slid::ParticleIDAlg::CalcPlaneLongitudinalDedx ( unsigned int  plane)
private

Calculate the longitudinal dE/dx in a given plane. Input is the actuall plane number in the detector.

Definition at line 506 of file ParticleIDAlg.cxx.

References CalcEnergyInLongitudinalPlane(), CalcTrkHitPath(), and fDedxLongByPlane.

Referenced by CalcPlaneLongDedxProb(), GetMIPPlanes(), and PlaneLongDedx().

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  }
double CalcEnergyInLongitudinalPlane(int iPlane)
Sum the energy in a longitudinal plane.
std::map< int, double > fDedxLongByPlane
double CalcTrkHitPath(unsigned int plane)
Calculate path length of shower through this plane.
double slid::ParticleIDAlg::CalcPlaneTransDedxProb ( DedxParticleType  particleType,
int  igev,
unsigned int  iTransversePlane 
)
private

Get transverse dE/dx probability for a given transverse plane, energy bin, and particle hypothesis.

Definition at line 713 of file ParticleIDAlg.cxx.

References CalcPlaneTransverseDedx(), fDedxDistributions, GetEntries(), HandyFuncs::Interpolate(), iReg, slid::kNumEnergyBin, slid::kNumTransversePlane, and slid::kNumXYRegion.

Referenced by CalcInterCellTransDedxProb().

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  }
static const int kNumTransversePlane
Number of "transverse" planes considered for dE/dx histograms ("planes" transverse to shower axis)...
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...
static const int kNumXYRegion
Number of XY regions into which detector is divided for dE/dx histograms.
std::vector< slid::DedxDistribution > fDedxDistributions
dE/dx distributions (histogram files and data) for each detector type and particle category ...
cout<< t1-> GetEntries()
Definition: plottest35.C:29
def Interpolate(x1, y1, x2, y2, yvalue)
Definition: HandyFuncs.py:218
static const int kNumEnergyBin
Number of energy bins into which detector is divided for dE/dx histograms.
double slid::ParticleIDAlg::CalcPlaneTransverseDedx ( unsigned int  iTransversePlane)
private

Calculate the transverse dE/dx in a given tranverse slice. Input is the transverse slice index.

Definition at line 530 of file ParticleIDAlg.cxx.

References CalcTrkHitPath(), slid::NuEEnergyAlg::CellEnergy(), dir, rb::Prong::Dir(), eventID, stan::math::fabs(), fDedxTransByPlane, fGeom, fNuEEnergyAlg, fPlaneCentroid, fPlaneHits, fShower, fShowerCol, fTransEnergy, fTransPathLength, makeTrainCVSamples::int, std::isnan(), geo::PlaneGeo::Ncells(), and geo::GeometryBase::Plane().

Referenced by CalcPlaneTransDedxProb(), DedxTransLL(), and PlaneTransDedx().

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  }
NuEEnergyAlg * fNuEEnergyAlg
Nue energy algorithms.
std::map< int, double > fTransEnergy
Map of transverse energy for each transverse plane.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
const PlaneGeo * Plane(unsigned int i) const
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
std::vector< const rb::Shower * > fShowerCol
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
art::ServiceHandle< geo::Geometry > fGeom
geometry
art::EventID eventID
data/MC flag
std::map< int, double > fTransPathLength
Map of path length for each transverse plane.
std::map< int, double > fDedxTransByPlane
TDirectory * dir
Definition: macro.C:5
std::map< int, double > fPlaneCentroid
std::map< int, std::map< int, int > > fPlaneHits
double CalcTrkHitPath(unsigned int plane)
Calculate path length of shower through this plane.
double slid::ParticleIDAlg::CalcTrkHitPath ( unsigned int  plane)
private

Calculate path length of shower through this plane.

Definition at line 473 of file ParticleIDAlg.cxx.

References geo::PlaneGeo::Cell(), rb::Prong::Dir(), e, stan::math::fabs(), fGeom, fPlaneHits, fShower, geo::CellGeo::HalfD(), Mag(), geo::GeometryBase::Plane(), util::pythag(), rb::Prong::Start(), and rb::Shower::Stop().

Referenced by CalcCellPlaneTransverseDedx(), CalcPlaneLongitudinalDedx(), and CalcPlaneTransverseDedx().

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  }
double HalfD() const
Definition: CellGeo.cxx:205
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
TVector3 Stop() const
Endpoint of the shower.
Definition: Shower.cxx:55
art::ServiceHandle< geo::Geometry > fGeom
geometry
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
float Mag() const
std::map< int, std::map< int, int > > fPlaneHits
Float_t e
Definition: plot.C:35
void slid::ParticleIDAlg::CalculateEnergyBinIndex ( )
private

Calculate energy bin indices and energy values from the histograms corresponding energy of shower to be used to interpolate between bins in energy.

Definition at line 440 of file ParticleIDAlg.cxx.

References e0, e1, fShowerEnergyGev, igev0, igev1, slid::kEnergyBins, and slid::kNumEnergyBin.

Referenced by SetShower().

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  }
double e0
Energy values of energy bins adjacent to shower. Calculated once and cached for a given shower/event ...
double fShowerEnergyGev
Shower energy.
static const int kNumEnergyBin
Number of energy bins into which detector is divided for dE/dx histograms.
static const float kEnergyBins[]
bin boundaries of energy
int igev0
Shower adjacent energy bins. Calculated once and cached for a given shower/event combo.
double slid::ParticleIDAlg::CellPlaneDedx ( int  tpIdx,
unsigned int  pIdx 
)

return transverse dedx for specified transverse plane

Definition at line 891 of file ParticleIDAlg.cxx.

References CalcCellPlaneTransverseDedx(), om::cout, allTimeWatchdog::endl, and slid::kNumTransversePlane.

Referenced by slid::LIDTraining::analyze(), and slid::LIDBuilder::produce().

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  }
static const int kNumTransversePlane
Number of "transverse" planes considered for dE/dx histograms ("planes" transverse to shower axis)...
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...
OStream cout
Definition: OStream.cxx:6
double slid::ParticleIDAlg::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 particle direction if opposite to the reconstructed direction.

Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis.

Definition at line 150 of file ParticleIDAlg.cxx.

References DedxLongLL(), rb::Prong::Dir(), rb::Prong::SetDir(), rb::Prong::SetStart(), and rb::Shower::Stop().

Referenced by slid::LIDBuilder::produce().

154  {
155  rb::Shower newShower = *vShower;
156  newShower.SetStart( vShower->Stop());
157  newShower.SetDir( -(vShower->Dir()) );
158  return DedxLongLL(partHyp, &(newShower), showercol, evt);
159  }
virtual void SetStart(TVector3 start)
Definition: Prong.cxx:75
virtual void SetDir(TVector3 dir)
Definition: Prong.cxx:104
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
TVector3 Stop() const
Endpoint of the shower.
Definition: Shower.cxx:55
A rb::Prong with a length.
Definition: Shower.h:18
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 slid::ParticleIDAlg::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.

Definition at line 101 of file ParticleIDAlg.cxx.

References CalcInterPlaneDedxProb(), rb::Prong::Dir(), fPartLongLL, fPlaneHits, test_ParserArtEvents::log, and SetShower().

Referenced by DedxInverseLongLL(), slid::LIDBuilder::produce(), and slid::SPIDBuilder::produce().

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  }
double CalcInterPlaneDedxProb(DedxParticleType type, unsigned int iLongPlane, unsigned int plane)
From Jianming&#39;s code, still working out changes.
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
std::map< int, std::map< int, int > > fPlaneHits
std::map< int, float > fPartLongLL
Map of the longitudinal ll by paricle type.
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.
double slid::ParticleIDAlg::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.

Definition at line 68 of file ParticleIDAlg.cxx.

References CalcInterCellTransDedxProb(), CalcPlaneTransverseDedx(), fPartTransLL, slid::kNumTransversePlane, test_ParserArtEvents::log, and SetShower().

Referenced by slid::LIDBuilder::produce(), and slid::SPIDBuilder::produce().

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  }
static const int kNumTransversePlane
Number of "transverse" planes considered for dE/dx histograms ("planes" transverse to shower axis)...
double CalcPlaneTransverseDedx(unsigned int iTransversePlane)
Calculate the transverse dE/dx in a given tranverse slice. Input is the transverse slice index...
std::map< int, float > fPartTransLL
Map of the transverse ll by particle type.
double CalcInterCellTransDedxProb(const DedxParticleType type, unsigned int iTransPlane)
From Jianming&#39;s code, still working out changes.
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.
TVector3 slid::ParticleIDAlg::GetCellNodePos ( unsigned int  plane1,
unsigned int  cell1,
unsigned int  plane2,
unsigned int  cell2 
)
private

Helper function for interia and asymmetry calculation.

Definition at line 1116 of file ParticleIDAlg.cxx.

References geo::PlaneGeo::Cell(), geom(), geo::CellGeo::GetCenter(), geo::kX, geo::kY, geo::GeometryBase::Plane(), geo::PlaneGeo::View(), submit_syst::x, submit_syst::y, and test::z.

Referenced by CalcAsymIneria().

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  }
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
Vertical planes which measure X.
Definition: PlaneGeo.h:28
const PlaneGeo * Plane(unsigned int i) const
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
z
Definition: test.py:28
void geom(int which=0)
Definition: geom.C:163
double slid::ParticleIDAlg::GetGapVertexToShowerStart ( const rb::Shower vShower,
TVector3  evtVtx,
const art::Event evt 
)

Calculate the gap between the vertex and the start of the shower.

Definition at line 767 of file ParticleIDAlg.cxx.

References Mag(), PointDoca(), std::sqrt(), rb::Prong::Start(), and rb::Shower::Stop().

Referenced by slid::LIDBuilder::produce().

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  }
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual TVector3 Start() const
Definition: Prong.h:73
TVector3 Stop() const
Endpoint of the shower.
Definition: Shower.cxx:55
double PointDoca(TVector3 vtx, TVector3 start, TVector3 stop)
Get distance of closest approach between shower and vertex.
float Mag() const
int slid::ParticleIDAlg::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.

Definition at line 1141 of file ParticleIDAlg.cxx.

References CalcPlaneLongitudinalDedx(), confusionMatrixTree::count, fMipRangeHigh, fMipRangeLow, fPlaneHits, and SetShower().

Referenced by slid::LIDBuilder::produce().

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
double fMipRangeLow
Lower bound of MIP Range.
double CalcPlaneLongitudinalDedx(unsigned int plane)
Calculate the longitudinal dE/dx in a given plane. Input is the actuall plane number in the detector...
std::map< int, std::map< int, int > > fPlaneHits
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.
TVector3 slid::ParticleIDAlg::InertiaXYZ ( ) const
inline

Return the inertia calculation.

Definition at line 134 of file ParticleIDAlg.h.

References febshutoff_auto::start.

134 {return fIxyz;}
TVector3 fIxyz
shower moment of ineria
bool slid::ParticleIDAlg::LoadDedxHistogramFiles ( const fhicl::ParameterSet pset,
novadaq::cnv::DetId  geom 
)
private

Load all of the histograms from files into arrays of TH1D* objects. The details of handling the histograms are taken care of within object of the helper class slid::DedxDistribution.

FD dE/dx input histogram files

NDOS dE/dx input histograms

ND dE/dx input histograms

Definition at line 306 of file ParticleIDAlg.cxx.

References om::cout, allTimeWatchdog::endl, util::EnvExpansion(), fDedxDistributions, fhicl::ParameterSet::get(), slid::kELECTRON, slid::kELECTRONCCCOH, slid::kELECTRONCCDIS, slid::kELECTRONCCQE, slid::kELECTRONCCRES, cmf::kFARDET, slid::kMUON, novadaq::cnv::kNDOS, cmf::kNEARDET, slid::kNEUTRON, slid::kPHOTON, slid::kPI0, slid::kPION, and slid::kPROTON.

Referenced by ParticleIDAlg().

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  }
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::vector< slid::DedxDistribution > fDedxDistributions
dE/dx distributions (histogram files and data) for each detector type and particle category ...
Prototype Near Detector on the surface at FNAL.
T get(std::string const &key) const
Definition: ParameterSet.h:231
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double slid::ParticleIDAlg::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.

Definition at line 799 of file ParticleIDAlg.cxx.

References e, evt, stan::math::fabs(), fClosestPi0Mass, fNuEEnergyAlg, fPi0Index, std::isnan(), calib::j, SetShower(), and slid::NuEEnergyAlg::ShowerEnergy().

Referenced by Pi0SecondPhoton(), and slid::LIDBuilder::produce().

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  }
NuEEnergyAlg * fNuEEnergyAlg
Nue energy algorithms.
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 isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
int evt
const double j
Definition: BetheBloch.cxx:29
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...
Float_t e
Definition: plot.C:35
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.
int slid::ParticleIDAlg::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 mass for the shower at index iShower.

Definition at line 835 of file ParticleIDAlg.cxx.

References fPi0Index, fShower, and Pi0Mass().

Referenced by slid::LIDBuilder::produce().

837  {
838  if( fShower == NULL)
839  Pi0Mass( showerCollection, iShower, evt);
840 
841  return fPi0Index;
842 
843  }// end of Pi0SecondPhoton
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
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...
double slid::ParticleIDAlg::PlaneCentroid ( unsigned int  pIdx)

Return the energy weighted mean shower positon for the plane in cm.

Definition at line 1024 of file ParticleIDAlg.cxx.

References fMapIdxToGlobalPlane, fPlaneCentroid, and it.

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  }
set< int >::iterator it
std::map< unsigned int, unsigned int > fMapIdxToGlobalPlane
std::map< int, double > fPlaneCentroid
double slid::ParticleIDAlg::PlaneHitCell ( unsigned int  pIdx)

Return point of intersection cell between shower core and plane.

Definition at line 968 of file ParticleIDAlg.cxx.

References geo::PlaneGeo::Cell(), om::cout, rb::Prong::Dir(), allTimeWatchdog::endl, stan::math::fabs(), fMapIdxToGlobalPlane, fShower, geom(), geo::CellGeo::GetCenter(), geo::CellGeo::HalfW(), geo::kX, geo::kY, geo::PlaneGeo::Ncells(), geo::GeometryBase::Plane(), NDAPDHVSetting::plane, rb::Prong::Start(), rb::Shower::Stop(), POTSpillRate::view, and geo::PlaneGeo::View().

Referenced by CalcCellPlaneTransverseDedx().

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  }
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
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
std::map< unsigned int, unsigned int > fMapIdxToGlobalPlane
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
TVector3 Stop() const
Endpoint of the shower.
Definition: Shower.cxx:55
OStream cout
Definition: OStream.cxx:6
void geom(int which=0)
Definition: geom.C:163
TVector3 slid::ParticleIDAlg::PlaneHitXYZ ( unsigned int  pIdx)

Return point of intersection between shower core and plane.

Definition at line 938 of file ParticleIDAlg.cxx.

References geo::PlaneGeo::Cell(), rb::Prong::Dir(), stan::math::fabs(), fMapIdxToGlobalPlane, fShower, geom(), geo::CellGeo::GetCenter(), geo::GeometryBase::Plane(), rb::Prong::Start(), and rb::Shower::Stop().

Referenced by slid::LIDTraining::analyze(), and CalcAsymIneria().

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  }
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
std::map< unsigned int, unsigned int > fMapIdxToGlobalPlane
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
TVector3 Stop() const
Endpoint of the shower.
Definition: Shower.cxx:55
void geom(int which=0)
Definition: geom.C:163
double slid::ParticleIDAlg::PlaneLongDedx ( unsigned int  pIdx)

return longitudinal dedx for a specified plane index

Definition at line 870 of file ParticleIDAlg.cxx.

References CalcPlaneLongitudinalDedx(), om::cout, allTimeWatchdog::endl, and fMapIdxToGlobalPlane.

Referenced by slid::LIDTraining::analyze(), slid::LIDBuilder::produce(), and slid::SPIDBuilder::produce().

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  }
std::map< unsigned int, unsigned int > fMapIdxToGlobalPlane
double CalcPlaneLongitudinalDedx(unsigned int plane)
Calculate the longitudinal dE/dx in a given plane. Input is the actuall plane number in the detector...
OStream cout
Definition: OStream.cxx:6
double slid::ParticleIDAlg::PlaneRadius ( unsigned int  pIdx)

return shower radius for a plane index

Definition at line 901 of file ParticleIDAlg.cxx.

References rb::CellHit::Cell(), rb::Cluster::Cell(), slid::NuEEnergyAlg::CellEnergy(), geo::ClosestApproach(), om::cout, dir, rb::Prong::Dir(), dist, allTimeWatchdog::endl, eventID, fGeom, fMapIdxToGlobalPlane, fNuEEnergyAlg, fPlaneHits, fShower, fShowerCol, it, geo::kX, rb::CellHit::Plane(), geo::GeometryBase::Plane(), Munits::rad, std::sqrt(), febshutoff_auto::start, rb::Prong::Start(), and rb::CellHit::View().

Referenced by slid::LIDTraining::analyze().

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  }
NuEEnergyAlg * fNuEEnergyAlg
Nue energy algorithms.
set< int >::iterator it
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
std::map< unsigned int, unsigned int > fMapIdxToGlobalPlane
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
double dist
Definition: runWimpSim.h:113
unsigned short Cell() const
Definition: CellHit.h:40
std::vector< const rb::Shower * > fShowerCol
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
art::ServiceHandle< geo::Geometry > fGeom
geometry
static constexpr Double_t rad
Definition: Munits.h:162
art::EventID eventID
data/MC flag
OStream cout
Definition: OStream.cxx:6
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 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
TDirectory * dir
Definition: macro.C:5
std::map< int, std::map< int, int > > fPlaneHits
unsigned int slid::ParticleIDAlg::PlaneToGlobalPlane ( unsigned int  pIdx)

Return the gloabel plane index for an index in the shower coordinate.

Definition at line 930 of file ParticleIDAlg.cxx.

References fMapIdxToGlobalPlane.

Referenced by slid::LIDTraining::analyze().

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  }
std::map< unsigned int, unsigned int > fMapIdxToGlobalPlane
double slid::ParticleIDAlg::PlaneTransDedx ( unsigned int  tpIdx)

return transverse dedx for specified transverse plane

Definition at line 881 of file ParticleIDAlg.cxx.

References CalcPlaneTransverseDedx(), om::cout, allTimeWatchdog::endl, and slid::kNumTransversePlane.

Referenced by slid::LIDTraining::analyze().

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  }
static const int kNumTransversePlane
Number of "transverse" planes considered for dE/dx histograms ("planes" transverse to shower axis)...
double CalcPlaneTransverseDedx(unsigned int iTransversePlane)
Calculate the transverse dE/dx in a given tranverse slice. Input is the transverse slice index...
OStream cout
Definition: OStream.cxx:6
double slid::ParticleIDAlg::PointDoca ( TVector3  vtx,
TVector3  start,
TVector3  stop 
)

Get distance of closest approach between shower and vertex.

Definition at line 786 of file ParticleIDAlg.cxx.

References d, Mag(), and febshutoff_auto::start.

Referenced by GetGapVertexToShowerStart(), slid::LIDBuilder::produce(), slid::SPIDBuilder::produce(), and cosrej::MakeNueCosRej::produce().

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  }
Float_t d
Definition: plot.C:236
float Mag() const
double slid::ParticleIDAlg::Radius ( const std::vector< const rb::Shower * >  showerCollection,
unsigned int  iShower,
const art::Event evt 
)

Calculate the shower radius.

Definition at line 846 of file ParticleIDAlg.cxx.

References rb::CellHit::Cell(), rb::Cluster::Cell(), slid::NuEEnergyAlg::CellEnergy(), geo::ClosestApproach(), dir, rb::Prong::Dir(), dist, eventID, evt, fGeom, fNuEEnergyAlg, fShower, fShowerCol, geo::kX, rb::Cluster::NCell(), rb::CellHit::Plane(), geo::GeometryBase::Plane(), Munits::rad, SetShower(), std::sqrt(), febshutoff_auto::start, rb::Prong::Start(), and rb::CellHit::View().

Referenced by slid::LIDBuilder::produce().

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  }
NuEEnergyAlg * fNuEEnergyAlg
Nue energy algorithms.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
const rb::Shower * fShower
rb::Shower object for which particle ID is to be calculated. This is cached.
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
double dist
Definition: runWimpSim.h:113
unsigned short Cell() const
Definition: CellHit.h:40
std::vector< const rb::Shower * > fShowerCol
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
int evt
art::ServiceHandle< geo::Geometry > fGeom
geometry
static constexpr Double_t rad
Definition: Munits.h:162
art::EventID eventID
data/MC flag
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 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
TDirectory * dir
Definition: macro.C:5
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.
void slid::ParticleIDAlg::reconfigure ( const fhicl::ParameterSet pset,
novadaq::cnv::DetId  detId 
)

Reinitialize parameters read from FCL file.

Definition at line 54 of file ParticleIDAlg.cxx.

56  {
57  }
void slid::ParticleIDAlg::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.

Set rb::Shower object to be analyzed. Calculate its total energy, which will be needed for calculations.

Parameters
vShower
evt

Definition at line 169 of file ParticleIDAlg.cxx.

References CalcAsymIneria(), CalcDetectorXYRegionIndex(), CalcPlaneCentroid(), CalcPlaneHits(), CalculateEnergyBinIndex(), e, gov::fnal::cd::rms::provider::equal(), eventID, fAsym, fDedxLongByPlane, fDedxTransByPlane, fIxyz, fMapIdxToGlobalPlane, fNuEEnergyAlg, fPartLongLL, fPartTransLL, fPi0Index, fPlaneCentroid, fPlaneHits, fShower, fShowerCol, fShowerEnergyGev, fTransEnergy, fTransPathLength, art::Event::id(), rb::Cluster::ID(), std::isnan(), and slid::NuEEnergyAlg::ShowerEnergy().

Referenced by slid::LIDTraining::analyze(), DedxLongLL(), DedxTransLL(), GetMIPPlanes(), Pi0Mass(), slid::LIDBuilder::produce(), slid::SPIDBuilder::produce(), and Radius().

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  }
NuEEnergyAlg * fNuEEnergyAlg
Nue energy algorithms.
std::map< int, double > fTransEnergy
Map of transverse energy for each transverse plane.
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...
bool equal(T *first, T *second)
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.
std::map< unsigned int, unsigned int > fMapIdxToGlobalPlane
double fShowerEnergyGev
Shower energy.
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
std::vector< const rb::Shower * > fShowerCol
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...
std::map< int, float > fPartTransLL
Map of the transverse ll by particle type.
void CalculateEnergyBinIndex()
Calculate energy bin indices and energy values from the histograms corresponding energy of shower to ...
double fAsym
shower asymmetry
art::EventID eventID
data/MC flag
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...
std::map< int, double > fTransPathLength
Map of path length for each transverse plane.
std::map< int, double > fDedxTransByPlane
TVector3 fIxyz
shower moment of ineria
const int ID() const
Definition: Cluster.h:75
std::map< int, double > fPlaneCentroid
std::map< int, std::map< int, int > > fPlaneHits
Float_t e
Definition: plot.C:35
std::map< int, float > fPartLongLL
Map of the longitudinal ll by paricle type.
EventID id() const
Definition: Event.h:56
std::map< int, double > fDedxLongByPlane
void CalcDetectorXYRegionIndex()
Determine detector XY region.
double slid::ParticleIDAlg::ShowerAsym ( ) const
inline

Return the asymmetry of the shower.

Definition at line 131 of file ParticleIDAlg.h.

131 {return fAsym;}
double fAsym
shower asymmetry

Member Data Documentation

double slid::ParticleIDAlg::e0
private

Energy values of energy bins adjacent to shower. Calculated once and cached for a given shower/event combo.

Definition at line 280 of file ParticleIDAlg.h.

Referenced by CalcInterCellTransDedxProb(), CalcInterPlaneDedxProb(), and CalculateEnergyBinIndex().

double slid::ParticleIDAlg::e1
private
art::EventID slid::ParticleIDAlg::eventID
private
double slid::ParticleIDAlg::fAsym
private

shower asymmetry

Definition at line 286 of file ParticleIDAlg.h.

Referenced by CalcAsymIneria(), and SetShower().

double slid::ParticleIDAlg::fClosestPi0Mass
private

Closest mass to pi0 mass formed by this shower and any other shower in the slice.

Definition at line 213 of file ParticleIDAlg.h.

Referenced by Pi0Mass().

std::vector<slid::DedxDistribution> slid::ParticleIDAlg::fDedxDistributions
private

dE/dx distributions (histogram files and data) for each detector type and particle category

Definition at line 216 of file ParticleIDAlg.h.

Referenced by CalcPlaneLongDedxProb(), CalcPlaneTransDedxProb(), and LoadDedxHistogramFiles().

std::map< int, double> slid::ParticleIDAlg::fDedxLongByPlane
private

Map of longitudinal dedx by plane. The first element is plane number and the second is the dedx for that plane. This is used for lazy caching of the dedx plane values in order to save recalculating for the same shower.

Definition at line 236 of file ParticleIDAlg.h.

Referenced by CalcPlaneLongitudinalDedx(), and SetShower().

std::map< int, double> slid::ParticleIDAlg::fDedxTransByPlane
private

Map of transverse dedx by plane. The first element is plane number and the second is the dedx for that plane. This is used for lazy caching of the dedx plane values in order to save recalculating for the same shower.

Definition at line 242 of file ParticleIDAlg.h.

Referenced by CalcPlaneTransverseDedx(), and SetShower().

novadaq::cnv::DetId slid::ParticleIDAlg::fDetId
private

Detector ID.

Definition at line 295 of file ParticleIDAlg.h.

Referenced by CalcDetectorXYRegionIndex().

art::ServiceHandle<geo::Geometry> slid::ParticleIDAlg::fGeom
private
TVector3 slid::ParticleIDAlg::fIxyz
private

shower moment of ineria

Definition at line 284 of file ParticleIDAlg.h.

Referenced by CalcAsymIneria(), and SetShower().

std::map<unsigned int, unsigned int> slid::ParticleIDAlg::fMapIdxToGlobalPlane
private

Map where the key is the index of the plane in the shower and the second is the global plane This map saves looping in some of the utility functions

Definition at line 253 of file ParticleIDAlg.h.

Referenced by CalcAsymIneria(), CalcCellPlaneTransverseDedx(), CalcPlaneHits(), PlaneCentroid(), PlaneHitCell(), PlaneHitXYZ(), PlaneLongDedx(), PlaneRadius(), PlaneToGlobalPlane(), and SetShower().

double slid::ParticleIDAlg::fMipRangeHigh
private

Upper bound of MIP Range.

Definition at line 289 of file ParticleIDAlg.h.

Referenced by GetMIPPlanes().

double slid::ParticleIDAlg::fMipRangeLow
private

Lower bound of MIP Range.

Definition at line 292 of file ParticleIDAlg.h.

Referenced by GetMIPPlanes().

NuEEnergyAlg* slid::ParticleIDAlg::fNuEEnergyAlg
private
std::map<int, float> slid::ParticleIDAlg::fPartLongLL

Map of the longitudinal ll by paricle type.

Definition at line 145 of file ParticleIDAlg.h.

Referenced by DedxLongLL(), slid::LIDBuilder::produce(), slid::SPIDBuilder::produce(), and SetShower().

std::map<int, float> slid::ParticleIDAlg::fPartTransLL

Map of the transverse ll by particle type.

Definition at line 148 of file ParticleIDAlg.h.

Referenced by DedxTransLL(), slid::LIDBuilder::produce(), slid::SPIDBuilder::produce(), and SetShower().

int slid::ParticleIDAlg::fPi0Index
private

Index of the other photon shower used for pi0 mass calculation for a given shower

Definition at line 267 of file ParticleIDAlg.h.

Referenced by Pi0Mass(), Pi0SecondPhoton(), and SetShower().

std::map< int, double > slid::ParticleIDAlg::fPlaneCentroid
private

Map where the first index is plane number and the second is the energy weighted center position for that plane

Definition at line 257 of file ParticleIDAlg.h.

Referenced by CalcCellPlaneTransverseDedx(), CalcPlaneTransverseDedx(), PlaneCentroid(), and SetShower().

std::map< int, std::map< int, int > > slid::ParticleIDAlg::fPlaneHits
private

Map where the first element is plane number and the second is a map with the key being the cell number in the plane and the second value is the global cell index in a shower

Definition at line 249 of file ParticleIDAlg.h.

Referenced by CalcAsymIneria(), CalcCellPlaneTransverseDedx(), CalcEnergyInLongitudinalPlane(), CalcPlaneCentroid(), CalcPlaneTransverseDedx(), CalcTrkHitPath(), DedxLongLL(), GetMIPPlanes(), PlaneRadius(), and SetShower().

const rb::Shower* slid::ParticleIDAlg::fShower
private
std::vector< const rb::Shower* > slid::ParticleIDAlg::fShowerCol
private
double slid::ParticleIDAlg::fShowerEnergyGev
private

Shower energy.

Definition at line 223 of file ParticleIDAlg.h.

Referenced by CalculateEnergyBinIndex(), and SetShower().

std::map<int, double > slid::ParticleIDAlg::fTransEnergy
private

Map of transverse energy for each transverse plane.

Definition at line 260 of file ParticleIDAlg.h.

Referenced by CalcPlaneTransverseDedx(), and SetShower().

std::map<int, double > slid::ParticleIDAlg::fTransPathLength
private

Map of path length for each transverse plane.

Definition at line 263 of file ParticleIDAlg.h.

Referenced by CalcPlaneTransverseDedx(), and SetShower().

int slid::ParticleIDAlg::igev0
private

Shower adjacent energy bins. Calculated once and cached for a given shower/event combo.

Definition at line 276 of file ParticleIDAlg.h.

Referenced by CalcInterCellTransDedxProb(), CalcInterPlaneDedxProb(), and CalculateEnergyBinIndex().

int slid::ParticleIDAlg::igev1
private
int slid::ParticleIDAlg::iReg
private

Detector sub-region bin. Calculated once and cached for a given shower/event combo.

Definition at line 270 of file ParticleIDAlg.h.

Referenced by CalcDetectorXYRegionIndex(), CalcPlaneLongDedxProb(), and CalcPlaneTransDedxProb().

TVector3 slid::ParticleIDAlg::pos
private

Shower start position.

Definition at line 273 of file ParticleIDAlg.h.

Referenced by CalcDetectorXYRegionIndex().


The documentation for this class was generated from the following files: