Public Types | Public Member Functions | Static Public Member Functions | Public Attributes | Protected Member Functions | Private Member Functions | Private Attributes | List of all members
jmshower::RecoJMShower Class Reference
Inheritance diagram for jmshower::RecoJMShower:
art::EDProducer art::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = ProducerBase::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 RecoJMShower (fhicl::ParameterSet const &pset)
 
 ~RecoJMShower ()
 
void preEvent (art::Event const &evt)
 
void postBeginRun (art::Run const &run)
 
void reconfigure (const fhicl::ParameterSet &pset)
 
void produce (art::Event &evt)
 
void beginRun (art::Run &run)
 
double GetPointDoca (TVector3 vtx, TVector3 start, TVector3 stop)
 
TVector3 GetCellDistToPoint (unsigned int plane, unsigned int cell, TVector3 point)
 
double GetCellDistToTrk (unsigned int plane, unsigned int cell, TVector3 start, TVector3 stop)
 
TVector3 GetShwStop (rb::Prong shw)
 
TVector3 GetShwStop (art::Ptr< rb::Prong > shw)
 
TVector3 GetTrkHitPos (unsigned int plane, unsigned int cell, art::Ptr< rb::Prong > trk)
 
TVector3 GetTrkHitPos (unsigned int plane, unsigned int cell, rb::Prong trk)
 
TVector3 GetTrkHitPos (unsigned int plane, unsigned int cell, TVector3 trkdir, TVector3 trkstart, TVector3 trkstop)
 
TVector3 GetTrkPlanePos (unsigned int plane, rb::Prong trk)
 
unsigned int GetTrkPlaneCell (unsigned int plane, art::Ptr< rb::Prong > trk)
 
unsigned int GetTrkPlaneCell (unsigned int plane, rb::Prong trk)
 
unsigned int GetTrkCPlaneCell (unsigned int plane, unsigned int i)
 
TVector3 GetTrkPlaneDistToEdge (unsigned int plane, art::Ptr< rb::Prong > trk)
 
TVector3 GetTrkPlaneDistToEdge (unsigned int plane, rb::Prong trk)
 
bool IsFiducial (rb::Prong shower)
 
double GetTrkDistToEdge (rb::Prong shower)
 
double GetTrkHitPath (unsigned int plane, unsigned int cell, rb::Prong shower)
 
unsigned int ContStartPlane (rb::Prong shower, unsigned int startp, int contp)
 
double DepositEnergy (rb::Cluster cluster, int i)
 
double DepositEnergy (rb::Prong shower)
 
double Energy (rb::Cluster cluster, int i)
 
double Energy (rb::Prong shower)
 
double ECorrMC (double gev)
 
double CellADCToGeV (double d, geo::View_t view, int detid)
 
double CellEDataMC (double gev, double w, geo::View_t view, int detid)
 
double PlaneDedxDataMC (double gev, int detid)
 
double CellTransDataMC (double gev, int detid)
 
double Radius (rb::Prong shower)
 
double Radius (rb::Cluster cluster, rb::Prong shower)
 
double RadiusV (rb::Cluster cluster)
 
double Gap (rb::Prong shower, TVector3 vertex)
 
double GetPlaneDedx (rb::Prong shower, unsigned int plane)
 
double GetPlaneDedxProb (rb::Prong shower, int ireg, int igev, unsigned int plane, double dedx, const int type)
 
double GetInterPlaneDedxProb (rb::Prong shower, unsigned int plane, double dedx, const int type, TVector3 pos, int iReg, int igev0, int igev1, double e0, double e1)
 
double GetDedxLongLL (rb::Prong shower, const int type)
 
double GetDedxInvLongLL (rb::Prong shower, const int type)
 
double GetDedxLongLL (rb::Prong shower, const int type, unsigned int startp, int contp)
 
double GetCellTransDedx (rb::Prong shower, unsigned int cell)
 
double GetCellTransDedxProb (rb::Prong shower, int ireg, int igev, unsigned int cell, double dedx, const int type)
 
double GetInterCellTransDedxProb (rb::Prong shower, unsigned int cell, double dedx, const int type, TVector3 pos, int iReg, int igev0, int igev1, double e0, double e1)
 
double GetDedxTransLL (rb::Prong shower, const int type)
 
bool IsMuon (jmshower::JMShower shower)
 
bool IsInvPhoton (jmshower::JMShower shower)
 
unsigned int NMIPPlane (rb::Prong shower)
 
unsigned int NMIPPlane (rb::Prong shower, unsigned int startp)
 
double CalcANN (unsigned int mode, jmshower::JMShower shower)
 
RVPInfo GetRVPStats (art::Handle< std::vector< rb::Cluster > > slices, unsigned int ic)
 
TVector3 GetCellNodePos (unsigned int plane1, unsigned int cell1, unsigned int plane2, unsigned int cell2)
 
TVector3 GetCentroid (rb::Cluster cluster)
 
int GetPlaneE1Cell (jmshower::PCluster pcluster)
 
int GetPlaneCentroidCell (jmshower::PCluster pcluster, rb::Prong shower)
 
std::vector< jmshower::JMShowerRecoShowers (art::EDProducer const &prod, art::Event &evt, std::unique_ptr< std::vector< jmshower::JMShower > > &recoshowercol, std::unique_ptr< art::Assns< jmshower::JMShower, rb::Prong > > &)
 
std::vector< jmshower::JMShowerExclShowers (const art::Event &evt)
 
template<typename PROD , BranchType B = InEvent>
ProductID getProductID (std::string const &instanceName={}) const
 
template<typename PROD , BranchType B>
ProductID getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
 
bool modifiesEvent () const
 
template<typename T , BranchType = InEvent>
ProductToken< Tconsumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< Tconsumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TconsumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< TmayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< TmayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TmayConsumeView (InputTag const &it)
 
base_engine_tcreateEngine (seed_t seed)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumernon_module_context ()
 

Public Attributes

TH1D * fHtPlaneDedxE [4][11][200]
 
TH1D * fHtPlaneDedxG [4][11][200]
 
TH1D * fHtPlaneDedxMu [4][11][200]
 
TH1D * fHtPlaneDedxPi0 [4][11][200]
 
TH1D * fHtPlaneDedxHad [4][11][200]
 
TH1D * fHtPlaneDedxP [4][11][200]
 
TH1D * fHtPlaneDedxN [4][11][200]
 
TH1D * fHtPlaneDedxPi [4][11][200]
 
TH1D * fHtPlaneDedxEqe [4][11][200]
 
TH1D * fHtPlaneDedxEres [4][11][200]
 
TH1D * fHtPlaneDedxEdis [4][11][200]
 
TH1D * fHtPlaneDedxEcoh [4][11][200]
 
TH1D * fHtPlaneDedxEsg [4][41][200]
 
unsigned fHtExpPlaneE [4][11]
 
unsigned fHtExpPlaneG [4][11]
 
unsigned fHtExpPlaneMu [4][11]
 
unsigned fHtExpPlanePi0 [4][11]
 
unsigned fHtExpPlaneHad [4][11]
 
unsigned fHtExpPlaneP [4][11]
 
unsigned fHtExpPlaneN [4][11]
 
unsigned fHtExpPlanePi [4][11]
 
unsigned fHtExpPlaneEqe [4][11]
 
unsigned fHtExpPlaneEres [4][11]
 
unsigned fHtExpPlaneEdis [4][11]
 
unsigned fHtExpPlaneEcoh [4][11]
 
unsigned fHtExpPlaneEsg [4][41]
 
TH1D * fHtCellTransDedxE [4][11][20]
 
TH1D * fHtCellTransDedxG [4][11][20]
 
TH1D * fHtCellTransDedxMu [4][11][20]
 
TH1D * fHtCellTransDedxPi0 [4][11][20]
 
TH1D * fHtCellTransDedxHad [4][11][20]
 
TH1D * fHtCellTransDedxP [4][11][20]
 
TH1D * fHtCellTransDedxN [4][11][20]
 
TH1D * fHtCellTransDedxPi [4][11][20]
 
TH1D * fHtCellTransDedxEqe [4][11][20]
 
TH1D * fHtCellTransDedxEres [4][11][20]
 
TH1D * fHtCellTransDedxEdis [4][11][20]
 
TH1D * fHtCellTransDedxEcoh [4][11][20]
 
TH1D * fHtCellTransDedxEsg [4][41][20]
 
std::string fFnameWeightsShape
 
std::string fFnameWeightsShapeElec
 
std::string fFnameWeightsShapeQE
 
std::string fFnameWeightsShapeRES
 
std::string fFnameWeightsShapeDIS
 
TMultiLayerPerceptron * mlp_shape
 
TMultiLayerPerceptron * mlp_shape_elec
 
TMultiLayerPerceptron * mlp_shape_E
 
TMultiLayerPerceptron * mlp_shape_NoCos
 
TMultiLayerPerceptron * mlp_shape_E_NoCos
 
TMultiLayerPerceptron * mlp_shape_QE
 
TMultiLayerPerceptron * mlp_shape_RES
 
TMultiLayerPerceptron * mlp_shape_DIS
 
TMultiLayerPerceptron * mlp_shape_EL
 
TMultiLayerPerceptron * mlp_shape_epi0
 
TMultiLayerPerceptron * mlp_shape_epi0EL
 
std::map< geo::OfflineChan, std::vector< double > > fHitEWeightMap
 
std::map< geo::OfflineChan, std::vector< double > > fHitDistMap
 
std::map< geo::OfflineChan, double > fHitEMap
 
std::vector< std::vector< jmshower::PCluster > > fShowerPClusterCol
 
std::vector< std::vector< jmshower::PCluster > > fTrackPClusterCol
 
std::vector< std::vector< jmshower::TCCluster > > fShowerTCClusterCol
 
std::vector< jmshower::JMShowerfExcludeShowerCol
 
std::string fLibPath
 
std::string pSliceDir
 
std::string pVertexLabel
 
std::string pInstLabel
 
std::string pTrackDir
 
double pCellHitPHThresh
 
int pPIDHist
 
int pRecoANN
 
bool pRecoPID
 
bool pRecoRVP
 
bool pRecoNode
 
bool pRecoPlaneR
 
bool pRecoInvLL
 
bool pRecluster
 
bool pRecoVtx
 
bool pVtxFit
 
bool pBadChannel
 
bool pCorAbsCell
 
bool pCorPigtail
 
bool pCorDataCell
 
bool pCorDataDedx
 
bool pUseUncalib
 
int pXEdgeId1
 
int pXEdgeId2
 
int pXEdgeId3
 
int pXEdgeId4
 
int pYEdgeId1
 
int pYEdgeId2
 
int pYEdgeId3
 
int pYEdgeId4
 
int pZEdgeId1
 
int pZEdgeId2
 
int pZEdgeId3
 
int pZEdgeId4
 
double pXEdge1
 
double pXEdge2
 
double pXEdge3
 
double pXEdge4
 
double pYEdge1
 
double pYEdge2
 
double pYEdge3
 
double pYEdge4
 
double pZEdge1
 
double pZEdge2
 
double pZEdge3
 
double pZEdge4
 

Protected Member Functions

CurrentProcessingContext const * currentContext () const
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

Private Member Functions

bool PreSel (const art::Event &evt)
 
bool LoadDedxHistograms (int detid)
 
bool LoadTreeWeights (int detid)
 

Private Attributes

int fDetid
 
bool isrealdata
 
bool fLoaded
 Have we done beginRun once yet? More...
 
std::vector< intfprongidx
 

Detailed Description

Definition at line 77 of file RecoJMShower_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 43 of file EDProducer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>
inherited

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

jmshower::RecoJMShower::RecoJMShower ( fhicl::ParameterSet const &  pset)

Definition at line 366 of file RecoJMShower_module.cc.

References reconfigure().

366  :
367  pSliceDir(pset.get<std::string>("SliceDir")),
368  pTrackDir(pset.get<std::string>("TrackDir")),
369  pCellHitPHThresh(50.),
370  pPIDHist(1),
371  pRecoANN(0),
372  pRecoPID(true),
373  pRecoRVP(true),
374  pRecoNode(true),
375  pRecoPlaneR(true),
376  pRecoInvLL(true),
377  pRecluster(true),
378  pRecoVtx(true),
379  pVtxFit(false),
380  pBadChannel(false),
381  pCorAbsCell(false),
382  pCorPigtail(true),
383  pCorDataCell(false),
384  pCorDataDedx(false),
385  pUseUncalib(false),
386  pXEdge1(-774.7),
387  pXEdge2(774.7),
388  pXEdge3(-774.7),
389  pXEdge4(-774.7),
390  pYEdge1(-774.7),
391  pYEdge2(774.7),
392  pYEdge3(-774.7),
393  pYEdge4(-774.7),
394  pZEdge1(0.0),
395  pZEdge2(5959.08),
396  pZEdge3(0.0),
397  pZEdge4(0.0),
398  fLoaded(false)
399  {
400  reconfigure(pset);
401  produces< std::vector<jmshower::JMShower> >();
402  produces< art::Assns<jmshower::JMShower, rb::Cluster> >();
403  produces< art::Assns<jmshower::JMShower, rb::Prong> >();
404  }
void reconfigure(const fhicl::ParameterSet &pset)
bool fLoaded
Have we done beginRun once yet?
enum BeamMode string
jmshower::RecoJMShower::~RecoJMShower ( )

Definition at line 407 of file RecoJMShower_module.cc.

408  {
409  }

Member Function Documentation

void jmshower::RecoJMShower::beginRun ( art::Run run)
virtual

Reimplemented from art::EDProducer.

Definition at line 464 of file RecoJMShower_module.cc.

References geo::GeometryBase::DetId(), fDetid, fLoaded, geom(), LoadDedxHistograms(), LoadTreeWeights(), and pRecoPID.

465  {
466  if(!fLoaded){
468  fDetid = geom->DetId();
469  if(pRecoPID==true)LoadDedxHistograms(geom->DetId());
470  LoadTreeWeights(geom->DetId());
471  fLoaded = true;
472  }
473 
474 
475  }
bool LoadDedxHistograms(int detid)
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
void geom(int which=0)
Definition: geom.C:163
bool fLoaded
Have we done beginRun once yet?
double jmshower::RecoJMShower::CalcANN ( unsigned int  mode,
jmshower::JMShower  shower 
)

Definition at line 5263 of file RecoJMShower_module.cc.

References std::cos(), om::cout, jmshower::JMShower::DedxLongLL(), jmshower::JMShower::DedxTransLL(), e, allTimeWatchdog::endl, jmshower::JMShower::Energy(), jmshower::JMShower::Gap(), GetPlaneDedx(), MECModelEnuComparisons::i, makeTrainCVSamples::int, jmshower::kElectron, jmshower::kElectronCOH, jmshower::kElectronDIS, jmshower::kElectronQE, jmshower::kElectronRES, jmshower::kElectronSG, jmshower::kGamma, jmshower::kMuon, jmshower::kNeutron, jmshower::kNull, jmshower::kPi0, jmshower::kPion, jmshower::kProton, std::max(), mlp_shape, mlp_shape_E, mlp_shape_E_NoCos, mlp_shape_EL, mlp_shape_elec, mlp_shape_epi0, mlp_shape_epi0EL, mlp_shape_NoCos, jmshower::JMShower::NPlane(), jmshower::JMShower::NueEnergy(), jmshower::JMShower::Pi0Mgg(), jmshower::JMShower::SliceEtot(), jmshower::JMShower::Theta(), and jmshower::JMShower::VtxGeV().

Referenced by RecoShowers().

5263  {
5264 // if(mode!=0&&mode!=1&&mode!=2&&&&mode!=3&&&&mode!=4) return -99;
5265  double ellmax = -9999;
5266  double elllmax = -9999;
5267  double elltmax = -9999;
5268  for(int itype = int(jmshower::kNull)+1; itype != int(jmshower::kElectronSG); ++itype){
5269  if(itype == int(jmshower::kElectron)||itype == int(jmshower::kElectronQE)||itype == int(jmshower::kElectronRES)||itype == int(jmshower::kElectronDIS)||itype == int(jmshower::kElectronCOH)){
5270  if(ellmax<shower.DedxLongLL(itype)+shower.DedxTransLL(itype)){
5271  ellmax = shower.DedxLongLL(itype)+shower.DedxTransLL(itype);
5272  elllmax = shower.DedxLongLL(itype);
5273  elltmax = shower.DedxTransLL(itype);
5274  }
5275  }
5276  }
5277 
5278  double egLLL = elllmax - shower.DedxLongLL(int(jmshower::kGamma));
5279  double egLLT = elltmax - shower.DedxTransLL(int(jmshower::kGamma));
5280  double emuLLL = elllmax - shower.DedxLongLL(int(jmshower::kMuon));
5281  double emuLLT = elltmax - shower.DedxTransLL(int(jmshower::kMuon));
5282  double epi0LLL = elllmax - shower.DedxLongLL(int(jmshower::kPi0));
5283  double epi0LLT = elltmax - shower.DedxTransLL(int(jmshower::kPi0));
5284  double epLLL = elllmax - shower.DedxLongLL(int(jmshower::kProton));
5285  double epLLT = elltmax - shower.DedxTransLL(int(jmshower::kProton));
5286  double enLLL = elllmax - shower.DedxLongLL(int(jmshower::kNeutron));
5287  double enLLT = elltmax - shower.DedxTransLL(int(jmshower::kNeutron));
5288  double epiLLL = elllmax - shower.DedxLongLL(int(jmshower::kPion));
5289  double epiLLT = elltmax - shower.DedxTransLL(int(jmshower::kPion));
5290 
5291  double eelgLLL = shower.DedxLongLL(int(jmshower::kElectronSG)) - shower.DedxLongLL(int(jmshower::kGamma));
5292  double eelgLLT = shower.DedxTransLL(int(jmshower::kElectronSG)) - shower.DedxTransLL(int(jmshower::kGamma));
5293  double eelmuLLL = shower.DedxLongLL(int(jmshower::kElectronSG)) - shower.DedxLongLL(int(jmshower::kMuon));
5294  double eelmuLLT = shower.DedxTransLL(int(jmshower::kElectronSG)) - shower.DedxTransLL(int(jmshower::kMuon));
5295  double eelpi0LLL = shower.DedxLongLL(int(jmshower::kElectronSG)) - shower.DedxLongLL(int(jmshower::kPi0));
5296  double eelpi0LLT = shower.DedxTransLL(int(jmshower::kElectronSG)) - shower.DedxTransLL(int(jmshower::kPi0));
5297  double eelpLLL = shower.DedxLongLL(int(jmshower::kElectronSG)) - shower.DedxLongLL(int(jmshower::kProton));
5298  double eelpLLT = shower.DedxTransLL(int(jmshower::kElectronSG)) - shower.DedxTransLL(int(jmshower::kProton));
5299  double eelnLLL = shower.DedxLongLL(int(jmshower::kElectronSG)) - shower.DedxLongLL(int(jmshower::kNeutron));
5300  double eelnLLT = shower.DedxTransLL(int(jmshower::kElectronSG)) - shower.DedxTransLL(int(jmshower::kNeutron));
5301  double eelpiLLL = shower.DedxLongLL(int(jmshower::kElectronSG)) - shower.DedxLongLL(int(jmshower::kPion));
5302  double eelpiLLT = shower.DedxTransLL(int(jmshower::kElectronSG)) - shower.DedxTransLL(int(jmshower::kPion));
5303 
5304  double vtxgev = shower.VtxGeV();
5305 // double vtxcentergev = shower.VtxCenterGeV();
5306 // double vtxcentergev5cell = shower.VtxCenterGeV5cell();
5307 // double vtxcentergev10cell = shower.VtxCenterGeV10cell();
5308 // double vtxcentergev15cell = shower.VtxCenterGeV15cell();
5309 // double vtxcentergev20cell = shower.VtxCenterGeV20cell();
5310 
5311  double pi0mass = std::max(shower.Pi0Mgg(),0.);
5312  double shE = 0;
5313  if(shower.SliceEtot()>1e-10)shE=shower.Energy()/shower.SliceEtot();//Use calibrated cell energy on showers.
5314 // if(shower.SliceGeV()>1e-10)shE=shower.Energy()/shower.SliceGeV();// Use total cell energy, calibrated by the slice mean position
5315  double nuE = shower.NueEnergy();
5316 // double nuE = (shower.Energy() + 0.282525 + 1.0766*(shower.SliceEtot()-shower.Energy()))/1.015;
5317 // double nuE = (shower.Energy() + 0.282525 + 1.0766*(shower.SliceGeV()-shower.DepositEnergy()))/1.015;
5318  double gap = shower.Gap();
5319  double costheta = cos(shower.Theta());
5320 // double length = 0;
5321 // if(shower.Energy()>1e-10)length = shower.TotalLength()/shower.Energy();
5322 
5323  Double_t params_shape[17];
5324  for(int i=0;i<17;i++){
5325  params_shape[i]=0;
5326  }
5327  params_shape[0] = egLLL;
5328  params_shape[1] = egLLT;
5329  params_shape[2] = emuLLL;
5330  params_shape[3] = emuLLT;
5331  params_shape[4] = epi0LLL;
5332  params_shape[5] = epi0LLT;
5333  params_shape[6] = epLLL;
5334  params_shape[7] = epLLT;
5335  params_shape[8] = enLLL;
5336  params_shape[9] = enLLT;
5337  params_shape[10] = epiLLL;
5338  params_shape[11] = epiLLT;
5339  params_shape[12] = pi0mass;
5340  params_shape[13] = shE;
5341  params_shape[14] = vtxgev;
5342  params_shape[15] = gap;
5343  params_shape[16] = costheta;
5344 
5345 
5346  Double_t params_shape_NoCos[16];
5347  for(int i=0;i<16;i++){
5348  params_shape_NoCos[i]=0;
5349  }
5350  params_shape_NoCos[0] = egLLL;
5351  params_shape_NoCos[1] = egLLT;
5352  params_shape_NoCos[2] = emuLLL;
5353  params_shape_NoCos[3] = emuLLT;
5354  params_shape_NoCos[4] = epi0LLL;
5355  params_shape_NoCos[5] = epi0LLT;
5356  params_shape_NoCos[6] = epLLL;
5357  params_shape_NoCos[7] = epLLT;
5358  params_shape_NoCos[8] = enLLL;
5359  params_shape_NoCos[9] = enLLT;
5360  params_shape_NoCos[10] = epiLLL;
5361  params_shape_NoCos[11] = epiLLT;
5362  params_shape_NoCos[12] = pi0mass;
5363  params_shape_NoCos[13] = shE;
5364  params_shape_NoCos[14] = vtxgev;
5365  params_shape_NoCos[15] = gap;
5366 
5367  Double_t params_elec[12];
5368  for(int i=0;i<12;i++){
5369  params_elec[i]=0;
5370  }
5371  params_elec[0] = egLLL;
5372  params_elec[1] = egLLT;
5373  params_elec[2] = epi0LLL;
5374  params_elec[3] = epi0LLT;
5375  params_elec[4] = epLLL;
5376  params_elec[5] = epLLT;
5377  params_elec[6] = enLLL;
5378  params_elec[7] = enLLT;
5379  params_elec[8] = epiLLL;
5380  params_elec[9] = epiLLT;
5381  params_elec[10] = shE;
5382  params_elec[11] = gap;
5383 
5384  Double_t params_shape_E[18];
5385  for(int i=0;i<18;i++){
5386  params_shape_E[i]=0;
5387  }
5388  params_shape_E[0] = egLLL;
5389  params_shape_E[1] = egLLT;
5390  params_shape_E[2] = emuLLL;
5391  params_shape_E[3] = emuLLT;
5392  params_shape_E[4] = epi0LLL;
5393  params_shape_E[5] = epi0LLT;
5394  params_shape_E[6] = epLLL;
5395  params_shape_E[7] = epLLT;
5396  params_shape_E[8] = enLLL;
5397  params_shape_E[9] = enLLT;
5398  params_shape_E[10] = epiLLL;
5399  params_shape_E[11] = epiLLT;
5400  params_shape_E[12] = pi0mass;
5401  params_shape_E[13] = shE;
5402  params_shape_E[14] = vtxgev;
5403  params_shape_E[15] = gap;
5404  params_shape_E[16] = costheta;
5405  params_shape_E[17] = nuE;
5406 
5407 
5408  Double_t params_shape_E_NoCos[17];
5409  for(int i=0;i<17;i++){
5410  params_shape_E_NoCos[i]=0;
5411  }
5412  params_shape_E_NoCos[0] = egLLL;
5413  params_shape_E_NoCos[1] = egLLT;
5414  params_shape_E_NoCos[2] = emuLLL;
5415  params_shape_E_NoCos[3] = emuLLT;
5416  params_shape_E_NoCos[4] = epi0LLL;
5417  params_shape_E_NoCos[5] = epi0LLT;
5418  params_shape_E_NoCos[6] = epLLL;
5419  params_shape_E_NoCos[7] = epLLT;
5420  params_shape_E_NoCos[8] = enLLL;
5421  params_shape_E_NoCos[9] = enLLT;
5422  params_shape_E_NoCos[10] = epiLLL;
5423  params_shape_E_NoCos[11] = epiLLT;
5424  params_shape_E_NoCos[12] = pi0mass;
5425  params_shape_E_NoCos[13] = shE;
5426  params_shape_E_NoCos[14] = vtxgev;
5427  params_shape_E_NoCos[15] = gap;
5428  params_shape_E_NoCos[16] = nuE;
5429 
5430  Double_t params_shape_EL[12];
5431  for(int i=0;i<12;i++){
5432  params_shape_EL[i]=0;
5433  }
5434  params_shape_EL[0] = eelgLLL;
5435  params_shape_EL[1] = eelgLLT;
5436  params_shape_EL[2] = eelmuLLL;
5437  params_shape_EL[3] = eelmuLLT;
5438  params_shape_EL[4] = eelpi0LLL;
5439  params_shape_EL[5] = eelpi0LLT;
5440  params_shape_EL[6] = eelpLLL;
5441  params_shape_EL[7] = eelpLLT;
5442  params_shape_EL[8] = eelnLLL;
5443  params_shape_EL[9] = eelnLLT;
5444  params_shape_EL[10] = eelpiLLL;
5445  params_shape_EL[11] = eelpiLLT;
5446 
5447  double dedx0 = 0;
5448  double dedx1 = 0;
5449  double dedx2 = 0;
5450  double dedx3 = 0;
5451 
5452  if(shower.NPlane()>0)dedx0 = RecoJMShower::GetPlaneDedx(shower, 0);
5453  if(shower.NPlane()>1)dedx1 = RecoJMShower::GetPlaneDedx(shower, 1);
5454  if(shower.NPlane()>2)dedx2 = RecoJMShower::GetPlaneDedx(shower, 2);
5455  if(shower.NPlane()>3)dedx3 = RecoJMShower::GetPlaneDedx(shower, 3);
5456 
5457  Double_t params_shape_epi0[5];
5458  for(int i=0;i<5;i++){
5459  params_shape_epi0[i]=0;
5460  }
5461  params_shape_epi0[0] = dedx0;
5462  params_shape_epi0[1] = dedx1;
5463  params_shape_epi0[2] = dedx2;
5464  params_shape_epi0[3] = dedx3;
5465  params_shape_epi0[4] = costheta;
5466 
5467  Double_t params_shape_epi0EL[4];
5468  for(int i=0;i<4;i++){
5469  params_shape_epi0EL[i]=0;
5470  }
5471  params_shape_epi0EL[0] = dedx0;
5472  params_shape_epi0EL[1] = dedx1;
5473  params_shape_epi0EL[2] = dedx2;
5474  params_shape_epi0EL[3] = dedx3;
5475 
5476 /*
5477 
5478  Double_t params_shape[18];
5479  for(int i=0;i<18;i++){
5480  params_shape[i]=0;
5481  }
5482  params_shape[0] = egLLL;
5483  params_shape[1] = egLLT;
5484  params_shape[2] = emuLLL;
5485  params_shape[3] = emuLLT;
5486  params_shape[4] = epi0LLL;
5487  params_shape[5] = epi0LLT;
5488  params_shape[6] = epLLL;
5489  params_shape[7] = epLLT;
5490  params_shape[8] = enLLL;
5491  params_shape[9] = enLLT;
5492  params_shape[10] = epiLLL;
5493  params_shape[11] = epiLLT;
5494  params_shape[12] = pi0mass;
5495  params_shape[13] = shE;
5496  params_shape[14] = gap;
5497  params_shape[15] = costheta;
5498  params_shape[16] = vtxgev;
5499  params_shape[17] = length;
5500 
5501 */
5502 
5503 /*
5504  Double_t params_shape[11];
5505  for(int i=0;i<11;i++){
5506  params_shape[i]=0;
5507  }
5508  params_shape[0] = egLLL;
5509  params_shape[1] = egLLT;
5510  params_shape[2] = epi0LLL;
5511  params_shape[3] = epi0LLT;
5512  params_shape[4] = epLLL;
5513  params_shape[5] = epLLT;
5514  params_shape[6] = enLLL;
5515  params_shape[7] = enLLT;
5516  params_shape[8] = epiLLL;
5517  params_shape[9] = epiLLT;
5518  params_shape[10] = pi0mass;
5519  params_shape[11] = vtxgev;
5520 */
5521  if(mode==0){
5522  double ann = mlp_shape->Evaluate(0,params_shape);
5523 // std::cout<<"mode, ann "<<mode<<" "<<ann<<std::endl;
5524  return ann;
5525  }
5526 
5527  if(mode==1){
5528  double ann = mlp_shape_elec->Evaluate(0,params_elec);
5529 // std::cout<<"mode, ann "<<mode<<" "<<ann<<std::endl;
5530  return ann;
5531  }
5532 
5533  if(mode==2){
5534  double ann = mlp_shape_E->Evaluate(0,params_shape_E);
5535 // std::cout<<"mode, ann "<<mode<<" "<<ann<<std::endl;
5536  return ann;
5537  }
5538 
5539  if(mode==3){
5540  double ann = mlp_shape_NoCos->Evaluate(0,params_shape_NoCos);
5541  return ann;
5542  }
5543 
5544  if(mode==4){
5545  double ann = mlp_shape_E_NoCos->Evaluate(0,params_shape_E_NoCos);
5546  return ann;
5547  }
5548 
5549  if(mode==5){
5550  double ann = mlp_shape_EL->Evaluate(0,params_shape_EL);
5551  return ann;
5552  }
5553 
5554  if(mode==6){
5555  double ann = mlp_shape_epi0->Evaluate(0,params_shape_epi0);
5556  std::cout<<"JM epi0 "<<params_shape_epi0[0]<<" "<<params_shape_epi0[1]<<" "<<params_shape_epi0[2]<<" "<<params_shape_epi0[3]<<" "<<params_shape_epi0[4]<<" "<<ann<<std::endl;
5557  return ann;
5558  }
5559 
5560  if(mode==7){
5561  double ann = mlp_shape_epi0EL->Evaluate(0,params_shape_epi0EL);
5562  return ann;
5563  }
5564 
5565  double ann0 = mlp_shape->Evaluate(0,params_shape);
5566  return ann0;
5567 /*
5568  double ann[3] = {-99,-99,-99};
5569  ann[0] = mlp_shape_QE->Evaluate(0,params_shape);
5570  ann[1] = mlp_shape_RES->Evaluate(0,params_shape);
5571  ann[2] = mlp_shape_DIS->Evaluate(0,params_shape);
5572 
5573  double annmax=-99;
5574  for(int ii=0;ii<3;ii++){
5575  if(annmax<ann[ii]){
5576  annmax = ann[ii];
5577  }
5578  }
5579 */
5580  }
T max(const caf::Proxy< T > &a, T b)
double NueEnergy() const
Definition: JMShower.cxx:263
TMultiLayerPerceptron * mlp_shape_elec
TMultiLayerPerceptron * mlp_shape_epi0EL
double Gap() const
Definition: JMShower.cxx:415
TMultiLayerPerceptron * mlp_shape
TMultiLayerPerceptron * mlp_shape_E
double GetPlaneDedx(rb::Prong shower, unsigned int plane)
unsigned int NPlane() const
Definition: JMShower.cxx:472
double DedxLongLL(int type) const
Definition: JMShower.cxx:672
TMultiLayerPerceptron * mlp_shape_NoCos
double Theta() const
Definition: JMShower.cxx:226
double VtxGeV() const
Definition: JMShower.cxx:156
double SliceEtot() const
Definition: JMShower.cxx:254
TMultiLayerPerceptron * mlp_shape_epi0
OStream cout
Definition: OStream.cxx:6
double DedxTransLL(int type) const
Definition: JMShower.cxx:708
double Pi0Mgg() const
Definition: JMShower.cxx:208
double Energy() const
Definition: JMShower.cxx:147
T cos(T number)
Definition: d0nt_math.hpp:78
Float_t e
Definition: plot.C:35
TMultiLayerPerceptron * mlp_shape_EL
TMultiLayerPerceptron * mlp_shape_E_NoCos
double jmshower::RecoJMShower::CellADCToGeV ( double  d,
geo::View_t  view,
int  detid 
)

Definition at line 3592 of file RecoJMShower_module.cc.

References e, novadaq::cnv::kFARDET, geo::kX, and geo::kY.

Referenced by RecoShowers().

3592  {
3593  double Att = 4.08050e+03;
3594  if(detid==novadaq::cnv::kFARDET){
3595  if(view==geo::kX){
3596  Att = 4.08574e+03+3.49097e+04*TMath::Exp(-2.67297e-03*d)
3597  -1.62442e+04*TMath::Exp(-3.14363e-02*d-(-3.93281e-05)*d*d-1.49793e-08*d*d*d);
3598  }else if(view == geo::kY){
3599  Att = 4.09347e+03+3.50384e+04*TMath::Exp(-2.68025e-03*d)
3600  -1.74280e+04*TMath::Exp(-3.08694e-02*d-(-3.86870e-05)*d*d-1.48737e-08*d*d*d);
3601  }
3602  }else{
3603  if(view==geo::kX){
3604  Att = 4.08574e+03+3.49097e+04*TMath::Exp(-2.67297e-03*d)
3605  -1.62442e+04*TMath::Exp(-3.14363e-02*d-(-3.93281e-05)*d*d-1.49793e-08*d*d*d);
3606  }else if(view == geo::kY){
3607  Att = 4.09347e+03+3.50384e+04*TMath::Exp(-2.68025e-03*d)
3608  -1.74280e+04*TMath::Exp(-3.08694e-02*d-(-3.86870e-05)*d*d-1.48737e-08*d*d*d);
3609  }
3610  }
3611  return Att;
3612  }
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Far Detector at Ash River, MN.
Float_t d
Definition: plot.C:236
Float_t e
Definition: plot.C:35
double jmshower::RecoJMShower::CellEDataMC ( double  gev,
double  w,
geo::View_t  view,
int  detid 
)

Definition at line 3615 of file RecoJMShower_module.cc.

References novadaq::cnv::kFARDET, geo::kX, geo::kY, and w.

Referenced by RecoShowers().

3615  {//
3616  double ecor = 1; //S12-11-16
3617  if(gev>0&&w>0){
3618  if(detid==novadaq::cnv::kFARDET){
3619  if(view == geo::kX){
3620  ecor = 1.60823+0.000519042*w-2.47044e-07*w*w;
3621  }else if(view == geo::kY){
3622  ecor = 1.3538+0.00051901*w-1.86311e-07*w*w;
3623  }
3624  }else{
3625  if(view == geo::kX){
3626  ecor = 1.60823+0.000519042*w-2.47044e-07*w*w;
3627  }else if(view == geo::kY){
3628  ecor = 1.3538+0.00051901*w-1.86311e-07*w*w;
3629  }
3630  }
3631  }
3632  if(ecor<0.0001)ecor = 1.0;
3633  return ecor;
3634  }
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Far Detector at Ash River, MN.
Float_t w
Definition: plot.C:20
double jmshower::RecoJMShower::CellTransDataMC ( double  gev,
int  detid 
)

Definition at line 3649 of file RecoJMShower_module.cc.

References novadaq::cnv::kFARDET.

Referenced by GetCellTransDedx().

3649  {//
3650  double ecor = 1;
3651  if(gev>0){
3652  if(detid==novadaq::cnv::kFARDET){
3653  ecor = 1.0/1.336;
3654  }else{
3655  ecor = 1.0/1.336;
3656  }
3657  }
3658  if(ecor<0.0001)ecor = 1.0;
3659  return ecor;
3660  }
Far Detector at Ash River, MN.
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::consumes ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::consumes ( InputTag const &  it)
inherited

Definition at line 146 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

147 {
148  if (!moduleContext_)
149  return ProductToken<T>::invalid();
150 
151  consumables_[BT].emplace_back(ConsumableType::Product,
152  TypeID{typeid(T)},
153  it.label(),
154  it.instance(),
155  it.process());
156  return ProductToken<T>{it};
157 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 161 of file Consumer.h.

References T.

162 {
163  if (!moduleContext_)
164  return;
165 
166  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
167 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 171 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

172 {
173  if (!moduleContext_)
174  return ViewToken<T>::invalid();
175 
176  consumables_[BT].emplace_back(ConsumableType::ViewElement,
177  TypeID{typeid(T)},
178  it.label(),
179  it.instance(),
180  it.process());
181  return ViewToken<T>{it};
182 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
unsigned int jmshower::RecoJMShower::ContStartPlane ( rb::Prong  shower,
unsigned int  startp,
int  contp 
)

Definition at line 3492 of file RecoJMShower_module.cc.

References fShowerPClusterCol, GetPlaneDedx(), MECModelEnuComparisons::i, rb::Cluster::ID(), NDAPDHVSetting::plane, and gen_flatrecord::size.

Referenced by RecoShowers().

3492  {
3493  int i = shower.ID();
3494  if(int(fShowerPClusterCol[i].size()-startp)<contp)return startp;
3495  unsigned int newstartp = startp;
3496  for(unsigned int plane= startp; plane<fShowerPClusterCol[i].size();plane++){
3497  int starttag = 1;
3498  for(int ic=0;ic<contp;ic++){
3499  if(plane+ic>fShowerPClusterCol[i].size()-1){
3500  starttag=0;
3501  break;
3502  }
3503 // double path = RecoJMShower::GetTrkHitPath(fShowerPClusterCol[i][plane+ic].Plane(), 0, shower);
3504 // double ep = RecoJMShower::Energy(fShowerPClusterCol[i][plane+ic],shower.ID());
3505 // double dedx = ep/path;
3506  double dedx = RecoJMShower::GetPlaneDedx(shower, plane+ic);
3507  if(dedx<0.0005){
3508  starttag=0;
3509  break;
3510  }
3511  }
3512  if(starttag==1){
3513  newstartp = plane;
3514  break;
3515  }
3516  }
3517  return newstartp;
3518  }
std::vector< std::vector< jmshower::PCluster > > fShowerPClusterCol
double GetPlaneDedx(rb::Prong shower, unsigned int plane)
const int ID() const
Definition: Cluster.h:75
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited
CurrentProcessingContext const* art::EDProducer::currentContext ( ) const
protectedinherited
double jmshower::RecoJMShower::DepositEnergy ( rb::Cluster  cluster,
int  i 
)

Definition at line 3526 of file RecoJMShower_module.cc.

References rb::CellHit::Cell(), rb::Cluster::Cell(), fHitEWeightMap, MECModelEnuComparisons::i, findDuplicateFiles::key, rb::Cluster::NCell(), and rb::CellHit::Plane().

Referenced by RecoShowers().

3526  {
3527  double gev = 0;
3528  for (unsigned int nh=0;nh<cluster.NCell();nh++) {
3529  art::Ptr<rb::CellHit> cellhit = cluster.Cell(nh);
3530  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
3531  double cellEnergy = fHitEWeightMap[key][i];
3532  gev += cellEnergy;
3533  }
3534  return gev;
3535  }
std::map< geo::OfflineChan, std::vector< double > > fHitEWeightMap
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
unsigned short Cell() const
Definition: CellHit.h:40
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
A (plane, cell) pair.
Definition: OfflineChan.h:17
double jmshower::RecoJMShower::DepositEnergy ( rb::Prong  shower)

Definition at line 3537 of file RecoJMShower_module.cc.

References rb::CellHit::Cell(), rb::Cluster::Cell(), fHitEWeightMap, MECModelEnuComparisons::i, rb::Cluster::ID(), findDuplicateFiles::key, rb::Cluster::NCell(), and rb::CellHit::Plane().

3537  {
3538  double gev = 0;
3539  int i = shower.ID();
3540  for (unsigned int nh=0;nh<shower.NCell();nh++) {
3541  art::Ptr<rb::CellHit> cellhit = shower.Cell(nh);
3542  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
3543  double cellEnergy = fHitEWeightMap[key][i];
3544  gev += cellEnergy;
3545  }
3546  return gev;
3547  }
std::map< geo::OfflineChan, std::vector< double > > fHitEWeightMap
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
unsigned short Cell() const
Definition: CellHit.h:40
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
A (plane, cell) pair.
Definition: OfflineChan.h:17
const int ID() const
Definition: Cluster.h:75
double jmshower::RecoJMShower::ECorrMC ( double  gev)

Definition at line 3584 of file RecoJMShower_module.cc.

References pCorAbsCell.

Referenced by Energy(), and RecoShowers().

3584  {//EM shower energy correction
3585 // double ecor = 5.63616e-01; //without pigtail
3586  double ecor = 5.68142e-01; //with pigtail
3587  if(pCorAbsCell==false)ecor = 5.62343e-01; //no mc cor
3588 // double ecor = 9.05891e-01*(5.53444e-01+9.92889e-02*gev-4.30034e-02*gev*gev-5.55264e-01*exp(-1.61264e+00*gev-1.07062e+01*gev*gev));
3589  return ecor;
3590  }
double jmshower::RecoJMShower::Energy ( rb::Cluster  cluster,
int  i 
)

Definition at line 3549 of file RecoJMShower_module.cc.

References rb::CellHit::Cell(), rb::Cluster::Cell(), ECorrMC(), fHitEWeightMap, MECModelEnuComparisons::i, findDuplicateFiles::key, rb::Cluster::NCell(), and rb::CellHit::Plane().

Referenced by GetDedxInvLongLL(), GetDedxLongLL(), GetDedxTransLL(), GetPlaneDedx(), and RecoShowers().

3549  {
3550  double gev = 0;
3551  for (unsigned int nh=0;nh<cluster.NCell();nh++) {
3552  art::Ptr<rb::CellHit> cellhit = cluster.Cell(nh);
3553  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
3554  double cellEnergy = fHitEWeightMap[key][i];
3555  gev += cellEnergy;
3556  }
3557  double ecor = RecoJMShower::ECorrMC(gev);
3558  gev = gev/ecor;
3559  return gev;
3560  }
std::map< geo::OfflineChan, std::vector< double > > fHitEWeightMap
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
unsigned short Cell() const
Definition: CellHit.h:40
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
A (plane, cell) pair.
Definition: OfflineChan.h:17
double jmshower::RecoJMShower::Energy ( rb::Prong  shower)

Definition at line 3562 of file RecoJMShower_module.cc.

References rb::CellHit::Cell(), rb::Cluster::Cell(), ECorrMC(), fHitEWeightMap, GetShwStop(), MECModelEnuComparisons::i, rb::Cluster::ID(), findDuplicateFiles::key, rb::Cluster::NCell(), pCorAbsCell, rb::CellHit::Plane(), and rb::Prong::Start().

3562  {
3563  double gev = 0;
3564  int i = shower.ID();
3565  for (unsigned int nh=0;nh<shower.NCell();nh++) {
3566  art::Ptr<rb::CellHit> cellhit = shower.Cell(nh);
3567  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
3568  double cellEnergy = fHitEWeightMap[key][i];
3569  gev += cellEnergy;
3570  }
3571  double ecor = RecoJMShower::ECorrMC(gev);
3572  double avgx = (shower.Start()[0]+RecoJMShower::GetShwStop(shower)[0])/2.;
3573  double avgy = (shower.Start()[1]+RecoJMShower::GetShwStop(shower)[1])/2.;
3574 // double poscor = (9.9899e-01+7.46983e-05*avgx+7.34449e-05*avgy);//S12-06-17
3575 // double poscor=1.+8.18868e-05*avgx+8.20582e-05*avgy;//without pigtail
3576  double poscor=1.+6.28813e-05*avgx+6.15165e-05*avgy; //with pigtail
3577  if(pCorAbsCell==false) poscor = 1.0+1.04866e-04*avgx+6.40668e-05*avgy;
3578  gev = gev/(ecor*poscor);
3579  //std::cout<<"JMShower Energy calc, depE: "<<showerDepE<<" ecor "<<ecor<<" avgx "<<avgx<<" avgy "<<avgy<<" poscor "<<poscor<<" showerE "<<gev<<std::endl;
3580 // gev = gev/ecor;
3581  return gev;
3582  }
std::map< geo::OfflineChan, std::vector< double > > fHitEWeightMap
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
TVector3 GetShwStop(rb::Prong shw)
unsigned short Plane() const
Definition: CellHit.h:39
virtual TVector3 Start() const
Definition: Prong.h:73
unsigned short Cell() const
Definition: CellHit.h:40
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
A (plane, cell) pair.
Definition: OfflineChan.h:17
const int ID() const
Definition: Cluster.h:75
std::vector< jmshower::JMShower > jmshower::RecoJMShower::ExclShowers ( const art::Event evt)

Definition at line 2882 of file RecoJMShower_module.cc.

References fExcludeShowerCol.

2882  {
2883  return fExcludeShowerCol;
2884  }
std::vector< jmshower::JMShower > fExcludeShowerCol
double jmshower::RecoJMShower::Gap ( rb::Prong  shower,
TVector3  vertex 
)

Definition at line 3754 of file RecoJMShower_module.cc.

References Mag(), slidt::showerStart, and rb::Prong::Start().

3754  {
3755  TVector3 showerStart=shower.Start();
3756  double gap = (vertex-showerStart).Mag();
3757  return gap;
3758  }
Definition: event.h:34
virtual TVector3 Start() const
Definition: Prong.h:73
float Mag() const
double showerStart[3]
seed_t art::EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited
TVector3 jmshower::RecoJMShower::GetCellDistToPoint ( unsigned int  plane,
unsigned int  cell,
TVector3  point 
)

Definition at line 2919 of file RecoJMShower_module.cc.

References std::abs(), geo::PlaneGeo::Cell(), geom(), geo::CellGeo::GetCenter(), geo::kX, geo::kY, geo::GeometryBase::Plane(), POTSpillRate::view, and geo::PlaneGeo::View().

Referenced by RadiusV().

2919  {
2920  Double_t cellXYZ[3]={0,0,0};
2922  geom->Plane(plane)->Cell(cell)->GetCenter(cellXYZ);
2923  geo::View_t view = geom->Plane(plane)->View();
2924  double celldistx = -99;
2925  double celldisty = -99;
2926  double celldistz = -99;
2927  if(view == geo::kX){
2928  celldistx = std::abs(cellXYZ[0]-point[0]);
2929  celldisty = 0;
2930  celldistz = std::abs(cellXYZ[2]-point[2]);
2931  }else if(view == geo::kY){
2932  celldistx = 0;
2933  celldisty = std::abs(cellXYZ[1]-point[1]);
2934  celldistz = std::abs(cellXYZ[2]-point[2]);
2935  }
2936  TVector3 celldist(celldistx,celldisty,celldistz);
2937  return celldist;
2938  }
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
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
void geom(int which=0)
Definition: geom.C:163
double jmshower::RecoJMShower::GetCellDistToTrk ( unsigned int  plane,
unsigned int  cell,
TVector3  start,
TVector3  stop 
)

Definition at line 2941 of file RecoJMShower_module.cc.

References std::abs(), geo::PlaneGeo::Cell(), geom(), geo::CellGeo::GetCenter(), geo::kX, geo::kY, RunSnowGlobes::mb, mc, fetch_tb_beamline_files::md, Munits::mg, plot_validation_datamc::p1, plot_validation_datamc::p2, geo::GeometryBase::Plane(), q2, std::sqrt(), POTSpillRate::view, geo::PlaneGeo::View(), x1, submit_syst::x2, y1, and submit_syst::y2.

Referenced by Radius(), and RecoShowers().

2941  {
2942  Double_t cellXYZ[3]={0,0,0};
2944  geom->Plane(plane)->Cell(cell)->GetCenter(cellXYZ);
2945  geo::View_t view = geom->Plane(plane)->View();
2946 
2947  double x1 = start[0];
2948  double y1 = start[1];
2949  double z1 = start[2];
2950  double p1 = (stop[0] - start[0])/(stop - start).Mag();
2951  double q1 = (stop[1] - start[1])/(stop - start).Mag();
2952  double r1 = (stop[2] - start[2])/(stop - start).Mag();
2953 
2954  double x2 = cellXYZ[0];
2955  double y2 = cellXYZ[1];
2956  double z2 = cellXYZ[2];
2957 
2958  double p2 = 1;
2959  double q2 = 0;
2960  double r2 = 0;
2961  if(view == geo::kX){
2962  p2 = 0;
2963  q2 = 1;
2964  r2 = 0;
2965  }else if(view == geo::kY){
2966  p2 = 1;
2967  q2 = 0;
2968  r2 = 0;
2969  }
2970  double ma = x2-x1;
2971  double mb = y2-y1;
2972  double mc = z2-z1;
2973  double md = p1;
2974  double me = q1;
2975  double mf = r1;
2976  double mg = p2;
2977  double mh = q2;
2978  double mi = r2;
2979 
2980  TMatrixD dDeNorm1(2,2);
2981  TMatrixD dDeNorm2(2,2);
2982  TMatrixD dDeNorm3(2,2);
2983  dDeNorm1[0][0] = p1;
2984  dDeNorm1[0][1] = q1;
2985  dDeNorm1[1][0] = p2;
2986  dDeNorm1[1][1] = q2;
2987  double det1 = dDeNorm1[0][0]*dDeNorm1[1][1] - dDeNorm1[0][1]*dDeNorm1[1][0];
2988  dDeNorm2[0][0] = q1;
2989  dDeNorm2[0][1] = r1;
2990  dDeNorm2[1][0] = q2;
2991  dDeNorm2[1][1] = r2;
2992  double det2 = dDeNorm2[0][0]*dDeNorm2[1][1] - dDeNorm2[0][1]*dDeNorm2[1][0];
2993  dDeNorm3[0][0] = r1;
2994  dDeNorm3[0][1] = p1;
2995  dDeNorm3[1][0] = r2;
2996  dDeNorm3[1][1] = p2;
2997  double det3 = dDeNorm3[0][0]*dDeNorm3[1][1] - dDeNorm3[0][1]*dDeNorm3[1][0];
2998  double det0 = ma*me*mi + mb*mf*mg + mc*md*mh - mc*me*mg - mb*md*mi - ma*mf*mh;
2999  double celldist = std::abs(det0)/sqrt(det1*det1+det2*det2+det3*det3);
3000  return celldist;
3001  }
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
Float_t y1[n_points_granero]
Definition: compare.C:5
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
Float_t x1[n_points_granero]
Definition: compare.C:5
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
const PlaneGeo * Plane(unsigned int i) const
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Double_t q2[12][num]
Definition: f2_nu.C:137
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
Definition: NueSkimmer.h:24
void geom(int which=0)
Definition: geom.C:163
static constexpr Double_t mg
Definition: Munits.h:210
TVector3 jmshower::RecoJMShower::GetCellNodePos ( unsigned int  plane1,
unsigned int  cell1,
unsigned int  plane2,
unsigned int  cell2 
)

Definition at line 5772 of file RecoJMShower_module.cc.

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 RecoShowers().

5772  {
5773  double cellXYZ1[3]={0,0,0};
5774  double cellXYZ2[3]={0,0,0};
5776  geom->Plane(plane1)->Cell(cell1)->GetCenter(cellXYZ1);
5777  geom->Plane(plane2)->Cell(cell2)->GetCenter(cellXYZ2);
5778  geo::View_t view1 = geom->Plane(plane1)->View();
5779  geo::View_t view2 = geom->Plane(plane2)->View();
5780  TVector3 nullpos(-9999,-9999,-9999);
5781  if(view1==view2)return nullpos;
5782  double x=0,y=0,z=0;
5783  if(view1 == geo::kX && view2 == geo::kY){
5784  x = cellXYZ1[0];
5785  y = cellXYZ2[1];
5786  }else if(view1 == geo::kY && view2 == geo::kX){
5787  x = cellXYZ2[0];
5788  y = cellXYZ1[1];
5789  }
5790  z = (cellXYZ1[2]+cellXYZ2[2])/2.;
5791  TVector3 nodepos(x,y,z);
5792  return nodepos;
5793  }
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 jmshower::RecoJMShower::GetCellTransDedx ( rb::Prong  shower,
unsigned int  cell 
)

Definition at line 4736 of file RecoJMShower_module.cc.

References std::abs(), rb::CellHit::Cell(), cellPos(), CellTransDataMC(), PandAna.Demos.tute_pid_validation::cid, dir, rb::Prong::Dir(), e, make_cached_def::e1, make_cached_def::e2, lem_server::ep, fDetid, fHitEWeightMap, fShowerPClusterCol, geom(), GetTrkHitPath(), MECModelEnuComparisons::i, rb::Cluster::ID(), makeTrainCVSamples::int, isrealdata, findDuplicateFiles::key, nc, geo::PlaneGeo::Ncells(), path, pCorDataDedx, rb::CellHit::Plane(), geo::GeometryBase::Plane(), NDAPDHVSetting::plane, and gen_flatrecord::size.

Referenced by GetDedxTransLL(), and RecoShowers().

4736  {// calculate transverse cell dE/dx based on centroid
4737 //if(cell!=0)return 0;
4738  int i = shower.ID();
4739  double totpath = 0 ;
4740  double ep = 0;
4741  for(unsigned int plane=0;plane<fShowerPClusterCol[i].size();plane++){
4742  unsigned int absplane = fShowerPClusterCol[i][plane].Plane();
4744  double sumPos = 0;
4745  double gev = 0;
4746  for (unsigned int nh=0;nh<fShowerPClusterCol[i][plane].NCell();nh++) {
4747  art::Ptr<rb::CellHit> cellhit = fShowerPClusterCol[i][plane].Cell(nh);
4748  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
4749  double cellEnergy = fHitEWeightMap[key][i];
4750  double cellPos = (double)cellhit->Cell() + 0.5;
4751  sumPos += cellEnergy*cellPos;
4752  gev += cellEnergy;
4753  }
4754  if(gev<0.00001)continue;
4755  double avgPos = sumPos/gev;
4756  int cid = (int)avgPos;// calculate centroid
4757  double wm = avgPos - (double)cid;//energy splitting for the centroid cell
4758  double wp = 1-wm; //energy splitting for the centroid cell
4759  if(cid<0)cid=0;
4760  if(cid>(int)geom->Plane(absplane)->Ncells()-1)cid=(int)geom->Plane(absplane)->Ncells()-1;
4761  double dir = 1.;
4762  if(std::abs(shower.Dir()[2])>0.000001)dir=std::abs(shower.Dir()[2]);
4763 // double dx=shower.Dir()[0];// Validate dir
4764 // double dy=shower.Dir()[1];
4765 // double dz=shower.Dir()[2];
4766 // if(geom->Plane(absplane)->View()==geo::kX&&dx*dx+dz*dz>1e-10)dir=std::abs(dz/sqrt(dx*dx+dz*dz));
4767 // if(geom->Plane(absplane)->View()==geo::kY&&dy*dy+dz*dz>1e-10)dir=std::abs(dz/sqrt(dy*dy+dz*dz));
4768  int ic = (int)cell + 1;
4769  double dph = ((double)ic/dir)+wm;//upper limit for cells included in a transeverse cell index
4770  double dpl = ((double)(ic-1)/dir)+wm;//lower limit for cells included in a transeverse cell index
4771  if(cid+(int)dph<(int)geom->Plane(plane)->Ncells()){
4772  double eps1 = dph-(int)dph;
4773  double eps2 = (int)dpl+1-dpl;
4774  geo::OfflineChan keye1(absplane, cid+(int)dph);
4775  geo::OfflineChan keye2(absplane, cid+(int)dpl);
4776  double e1 = 0;
4777  double e2 = 0;
4778  double e = 0;
4779  if((int)fHitEWeightMap[keye1].size()>i)e1 = eps1*fHitEWeightMap[keye1][i];
4780  if((int)fHitEWeightMap[keye2].size()>i)e2 = eps2*fHitEWeightMap[keye2][i];
4781  int nc = 0;
4782  if((int)dph-(int)dpl>1){
4783  for(int ii=(int)dph-1;ii>=(int)dpl+1;ii--){
4784  geo::OfflineChan keye(absplane,cid+ii);
4785  if((int)fHitEWeightMap[keye].size()>i)e+=fHitEWeightMap[keye][i];
4786  nc++;
4787  }
4788  }
4789  double path = RecoJMShower::GetTrkHitPath(fShowerPClusterCol[i][plane].Plane(), 0, shower);
4790  totpath+=path;
4791  ep+=e1+e2+e;
4792 // Simple calculation by using cellE/path
4793 // double ecid = 0;
4794 // geo::OfflineChan keycid(absplane,cid);
4795 // if((int)fHitEWeightMap[keycid].size()>i)ecid = fHitEWeightMap[keycid][i];
4796 // double dir3 = 1;
4797 // if(std::abs(shower.Dir()[2])>0.000001)dir3=shower.Dir()[2];
4798 // ep+=ecid/std::abs(dir3)/2.;
4799  }
4800  double dmh = ((double)ic/dir)+wp;//upper limit for cells included in a transeverse cell index
4801  double dml = ((double)(ic-1)/dir)+wp;//lower limit for cells included in a transeverse cell index
4802 
4803  if(cid-(int)dmh>=0){
4804  double eps1 = dmh-(int)dmh;
4805  double eps2 = (int)dml+1-dml;
4806  geo::OfflineChan keye1(absplane, cid-(int)dmh);
4807  geo::OfflineChan keye2(absplane, cid-(int)dml);
4808  double e1 = 0;
4809  double e2 = 0;
4810  double e = 0;
4811  if((int)fHitEWeightMap[keye1].size()>i)e1 = eps1*fHitEWeightMap[keye1][i];
4812  if((int)fHitEWeightMap[keye2].size()>i)e2 = eps2*fHitEWeightMap[keye2][i];
4813  int nc = 0;
4814  if((int)dmh-(int)dml>1){
4815  for(int ii=(int)dmh-1;ii>=(int)dml+1;ii--){
4816  geo::OfflineChan keye(absplane,cid-ii);
4817  if((int)fHitEWeightMap[keye].size()>i)e+=fHitEWeightMap[keye][i];
4818  nc++;
4819  }
4820  }
4821  ep+=e1+e2+e;
4822 
4823 // Simple calculation by using cellE/path
4824 // double ecid = 0;
4825 // geo::OfflineChan keycid(absplane,cid);
4826 // if(fHitEWeightMap[keycid].size()>i)ecid = fHitEWeightMap[keycid][i];
4827 // double dir3 = 1;
4828 // if(std::abs(shower.Dir()[2])>0.000001)dir3=shower.Dir()[2];
4829 // ep+=ecid/std::abs(dir3)/2.;
4830 
4831  double path = RecoJMShower::GetTrkHitPath(fShowerPClusterCol[i][plane].Plane(), 0, shower);
4832  totpath+=path;
4833  }
4834  }
4835 // double ep = RecoJMShower::DepositEnergy(fShowerTCClusterCol[i][cell], i);
4836  double dedx = -1;
4837  if(totpath>0.00001)dedx = ep/(totpath/2.);
4838  if(isrealdata==true&&pCorDataDedx)dedx=dedx/RecoJMShower::CellTransDataMC(dedx, fDetid);
4839  //std::cout<<"jm trans "<<cell<<" p "<<totpath/2.<<" ep "<<ep<<" dedx "<<dedx<<std::endl;
4840  return dedx;
4841  }
std::map< geo::OfflineChan, std::vector< double > > fHitEWeightMap
unsigned short Plane() const
Definition: CellHit.h:39
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
const PlaneGeo * Plane(unsigned int i) const
std::vector< std::vector< jmshower::PCluster > > fShowerPClusterCol
float abs(float number)
Definition: d0nt_math.hpp:39
unsigned short Cell() const
Definition: CellHit.h:40
double CellTransDataMC(double gev, int detid)
double GetTrkHitPath(unsigned int plane, unsigned int cell, rb::Prong shower)
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
void cellPos()
Definition: cellShifts.C:26
const std::string path
Definition: plot_BEN.C:43
TDirectory * dir
Definition: macro.C:5
A (plane, cell) pair.
Definition: OfflineChan.h:17
void geom(int which=0)
Definition: geom.C:163
enum BeamMode nc
const int ID() const
Definition: Cluster.h:75
Float_t e
Definition: plot.C:35
double jmshower::RecoJMShower::GetCellTransDedxProb ( rb::Prong  shower,
int  ireg,
int  igev,
unsigned int  cell,
double  dedx,
const int  type 
)

Definition at line 4843 of file RecoJMShower_module.cc.

References getBrightness::cell, fHtCellTransDedxE, fHtCellTransDedxEcoh, fHtCellTransDedxEdis, fHtCellTransDedxEqe, fHtCellTransDedxEres, fHtCellTransDedxEsg, fHtCellTransDedxG, fHtCellTransDedxMu, fHtCellTransDedxN, fHtCellTransDedxP, fHtCellTransDedxPi, fHtCellTransDedxPi0, GetEntries(), HandyFuncs::Interpolate(), jmshower::kElectron, jmshower::kElectronCOH, jmshower::kElectronDIS, jmshower::kElectronQE, jmshower::kElectronRES, jmshower::kElectronSG, jmshower::kGamma, jmshower::kMuon, jmshower::kNeutron, jmshower::kPi0, jmshower::kPion, jmshower::kProton, and getGoodRuns4SAM::n.

Referenced by GetInterCellTransDedxProb().

4843  {
4844  if(int(jmshower::kElectronSG)==type){
4845  if(ireg<0||ireg>3)return 1;
4846  if(igev>=41||igev<0)return 1;
4847  if(cell>=20)return 1;
4848  }else{
4849  if(ireg<0||ireg>3)return 1;
4850  if(igev>=11||igev<0)return 1;
4851  if(cell>=20)return 1;
4852  }
4853 
4854 // if(dedx<=0.00003)return 1;
4855  double p = 0;
4856  if(int(jmshower::kElectron)==type){
4857  double ntot = double(fHtCellTransDedxE[ireg][igev][cell]->GetEntries());
4858  double n = double(fHtCellTransDedxE[ireg][igev][cell]->Interpolate(dedx));
4859 // if(dedx<=0.00003)n = double(fHtCellTransDedxE[ireg][igev][cell]->GetBinContent(1));
4860  double nbin = fHtCellTransDedxE[ireg][igev][cell]->GetNbinsX();
4861  if(n<0.0001)n=0.0001;
4862  if(ntot<1){
4863  p = 0;
4864  }else{
4865  p = n*nbin/ntot;
4866  }
4867  }
4868  if(int(jmshower::kGamma)==type){
4869  double ntot = double(fHtCellTransDedxG[ireg][igev][cell]->GetEntries());
4870  double n = double(fHtCellTransDedxG[ireg][igev][cell]->Interpolate(dedx));
4871 // if(dedx<=0.00003)n = double(fHtCellTransDedxG[ireg][igev][cell]->GetBinContent(1));
4872  double nbin = fHtCellTransDedxG[ireg][igev][cell]->GetNbinsX();
4873  if(n<0.0001)n=0.0001;
4874  if(ntot<1){
4875  p = 0;
4876  }else{
4877  p = n*nbin/ntot;
4878  }
4879  }
4880 
4881  if(int(jmshower::kMuon)==type){
4882  double ntot = double(fHtCellTransDedxMu[ireg][igev][cell]->GetEntries());
4883  double n = double(fHtCellTransDedxMu[ireg][igev][cell]->Interpolate(dedx));
4884  // if(dedx<=0.00003)n = double(fHtCellTransDedxMu[ireg][igev][cell]->GetBinContent(1));
4885  double nbin = fHtCellTransDedxMu[ireg][igev][cell]->GetNbinsX();
4886  if(n<0.0001)n=0.0001;
4887  if(ntot<1){
4888  p = 0;
4889  }else{
4890  p = n*nbin/ntot;
4891  }
4892  }
4893 
4894  if(int(jmshower::kPi0)==type){
4895  double ntot = double(fHtCellTransDedxPi0[ireg][igev][cell]->GetEntries());
4896  double n = double(fHtCellTransDedxPi0[ireg][igev][cell]->Interpolate(dedx));
4897 // if(dedx<=0.00003)n = double(fHtCellTransDedxPi0[ireg][igev][cell]->GetBinContent(1));
4898 // double n = double(fHtCellTransDedxPi0[ireg][igev][cell]->GetBinContent(fHtCellTransDedxPi0[ireg][igev][cell]->FindBin(dedx)));
4899 // int nnbin = fHtCellTransDedxPi0[ireg][igev][cell]->FindBin(dedx);
4900  double nbin = fHtCellTransDedxPi0[ireg][igev][cell]->GetNbinsX();
4901  if(n<0.0001)n=0.0001;
4902  if(ntot<1){
4903  p = 0;
4904  }else{
4905  p = n*nbin/ntot;
4906  }
4907  }
4908 /*
4909  if(int(jmshower::kHad)==type){
4910  double ntot = double(fHtCellTransDedxHad[ireg][igev][cell]->GetEntries());
4911  double n = double(fHtCellTransDedxHad[ireg][igev][cell]->Interpolate(dedx));
4912 // if(dedx<=0.00003)n = double(fHtCellTransDedxHad[ireg][igev][cell]->GetBinContent(1));
4913 // double n = double(fHtCellTransDedxHad[ireg][igev][cell]->GetBinContent(fHtCellTransDedxHad[ireg][igev][cell]->FindBin(dedx)));
4914 // int nnbin = fHtCellTransDedxHad[ireg][igev][cell]->FindBin(dedx);
4915  double nbin = fHtCellTransDedxHad[ireg][igev][cell]->GetNbinsX();
4916  if(n<0.0001)n=0.0001;
4917  if(ntot<1){
4918  p = 0;
4919  }else{
4920  p = n*nbin/ntot;
4921  }
4922  }
4923 */
4924 
4925 
4926 
4927  if(int(jmshower::kProton)==type){
4928  double ntot = double(fHtCellTransDedxP[ireg][igev][cell]->GetEntries());
4929  double n = double(fHtCellTransDedxP[ireg][igev][cell]->Interpolate(dedx));
4930 // if(dedx<=0.00003)n = double(fHtCellTransDedxP[ireg][igev][cell]->GetBinContent(1));
4931 // double n = double(fHtCellTransDedxP[ireg][igev][cell]->GetBinContent(fHtCellTransDedxP[ireg][igev][cell]->FindBin(dedx)));
4932 // int nnbin = fHtCellTransDedxP[ireg][igev][cell]->FindBin(dedx);
4933  double nbin = fHtCellTransDedxP[ireg][igev][cell]->GetNbinsX();
4934  if(n<0.0001)n=0.0001;
4935  if(ntot<1){
4936  p = 0;
4937  }else{
4938  p = n*nbin/ntot;
4939  }
4940  }
4941 
4942 
4943  if(int(jmshower::kNeutron)==type){
4944  double ntot = double(fHtCellTransDedxN[ireg][igev][cell]->GetEntries());
4945  double n = double(fHtCellTransDedxN[ireg][igev][cell]->Interpolate(dedx));
4946 // if(dedx<=0.00003)n = double(fHtCellTransDedxN[ireg][igev][cell]->GetBinContent(1));
4947 // double n = double(fHtCellTransDedxN[ireg][igev][cell]->GetBinContent(fHtCellTransDedxN[ireg][igev][cell]->FindBin(dedx)));
4948 // int nnbin = fHtCellTransDedxN[ireg][igev][cell]->FindBin(dedx);
4949  double nbin = fHtCellTransDedxN[ireg][igev][cell]->GetNbinsX();
4950  if(n<0.0001)n=0.0001;
4951  if(ntot<1){
4952  p = 0;
4953  }else{
4954  p = n*nbin/ntot;
4955  }
4956  }
4957 
4958  if(int(jmshower::kPion)==type){
4959  double ntot = double(fHtCellTransDedxPi[ireg][igev][cell]->GetEntries());
4960  double n = double(fHtCellTransDedxPi[ireg][igev][cell]->Interpolate(dedx));
4961 // if(dedx<=0.00003)n = double(fHtCellTransDedxPi[ireg][igev][cell]->GetBinContent(1));
4962 // double n = double(fHtCellTransDedxPi[ireg][igev][cell]->GetBinContent(fHtCellTransDedxPi[ireg][igev][cell]->FindBin(dedx)));
4963 // int nnbin = fHtCellTransDedxPi[ireg][igev][cell]->FindBin(dedx);
4964  double nbin = fHtCellTransDedxPi[ireg][igev][cell]->GetNbinsX();
4965  if(n<0.0001)n=0.0001;
4966  if(ntot<1){
4967  p = 0;
4968  }else{
4969  p = n*nbin/ntot;
4970  }
4971  }
4972 
4973  if(int(jmshower::kElectronQE)==type){
4974  double ntot = double(fHtCellTransDedxEqe[ireg][igev][cell]->GetEntries());
4975  double n = double(fHtCellTransDedxEqe[ireg][igev][cell]->Interpolate(dedx));
4976  double nbin = fHtCellTransDedxEqe[ireg][igev][cell]->GetNbinsX();
4977  if(n<0.0001)n=0.0001;
4978  if(ntot<1){
4979  p = 0;
4980 
4981  }else{
4982  p = n*nbin/ntot;
4983  }
4984  }
4985 
4986  if(int(jmshower::kElectronRES)==type){
4987  double ntot = double(fHtCellTransDedxEres[ireg][igev][cell]->GetEntries());
4988  double n = double(fHtCellTransDedxEres[ireg][igev][cell]->Interpolate(dedx));
4989  double nbin = fHtCellTransDedxEres[ireg][igev][cell]->GetNbinsX();
4990  if(n<0.0001)n=0.0001;
4991  if(ntot<1){
4992  p = 0;
4993  }else{
4994  p = n*nbin/ntot;
4995  }
4996  }
4997 
4998  if(int(jmshower::kElectronDIS)==type){
4999  double ntot = double(fHtCellTransDedxEdis[ireg][igev][cell]->GetEntries());
5000  double n = double(fHtCellTransDedxEdis[ireg][igev][cell]->Interpolate(dedx));
5001  double nbin = fHtCellTransDedxEdis[ireg][igev][cell]->GetNbinsX();
5002  if(n<0.0001)n=0.0001;
5003  if(ntot<1){
5004  p = 0;
5005  }else{
5006  p = n*nbin/ntot;
5007  }
5008  }
5009 
5010  if(int(jmshower::kElectronCOH)==type){
5011  double ntot = double(fHtCellTransDedxEcoh[ireg][igev][cell]->GetEntries());
5012  double n = double(fHtCellTransDedxEcoh[ireg][igev][cell]->Interpolate(dedx));
5013  double nbin = fHtCellTransDedxEcoh[ireg][igev][cell]->GetNbinsX();
5014  if(n<0.0001)n=0.0001;
5015  if(ntot<1){
5016  p = 0;
5017  }else{
5018  p = n*nbin/ntot;
5019  }
5020  }
5021 
5022  if(int(jmshower::kElectronSG)==type){
5023  double ntot = double(fHtCellTransDedxEsg[ireg][igev][cell]->GetEntries());
5024  double n = double(fHtCellTransDedxEsg[ireg][igev][cell]->Interpolate(dedx));
5025  double nbin = fHtCellTransDedxEsg[ireg][igev][cell]->GetNbinsX();
5026  if(n<0.0001)n=0.0001;
5027  if(ntot<1){
5028  p = 0;
5029  }else{
5030  p = n*nbin/ntot;
5031  }
5032  }
5033 
5034  return p;
5035  }
TH1D * fHtCellTransDedxPi[4][11][20]
const char * p
Definition: xmltok.h:285
cout<< t1-> GetEntries()
Definition: plottest35.C:29
TH1D * fHtCellTransDedxE[4][11][20]
TH1D * fHtCellTransDedxP[4][11][20]
TH1D * fHtCellTransDedxG[4][11][20]
TH1D * fHtCellTransDedxPi0[4][11][20]
TH1D * fHtCellTransDedxN[4][11][20]
def Interpolate(x1, y1, x2, y2, yvalue)
Definition: HandyFuncs.py:218
TH1D * fHtCellTransDedxMu[4][11][20]
TH1D * fHtCellTransDedxEres[4][11][20]
TH1D * fHtCellTransDedxEdis[4][11][20]
TH1D * fHtCellTransDedxEcoh[4][11][20]
TH1D * fHtCellTransDedxEqe[4][11][20]
TH1D * fHtCellTransDedxEsg[4][41][20]
TVector3 jmshower::RecoJMShower::GetCentroid ( rb::Cluster  cluster)

Definition at line 3712 of file RecoJMShower_module.cc.

References rawdata::RawDigit::ADC(), rb::CellHit::Cell(), geom(), rb::Cluster::NXCell(), rb::Cluster::NYCell(), rb::CellHit::Plane(), geo::GeometryBase::Plane(), rb::Cluster::XCell(), and rb::Cluster::YCell().

Referenced by RadiusV(), and RecoShowers().

3712  {
3713  double xadc = 0;
3714  double yadc = 0;
3715  double zadc = 0;
3716  double xadccellx = 0;
3717  double yadccelly = 0;
3718  double zadccellz = 0;
3719  double cx = 0;
3720  double cy = 0;
3721  double cz = 0;
3722 
3723  for(unsigned int nh=0; nh<cluster.NXCell();nh++){
3724  art::Ptr<rb::CellHit> cellhit = cluster.XCell(nh);
3725  Double_t cellXYZ[3]={0,0,0};
3727  geom->Plane(cellhit->Plane())->Cell(cellhit->Cell())->GetCenter(cellXYZ);
3728  double cellADC = double(cellhit->ADC());
3729  xadc += cellADC;
3730  zadc += cellADC;
3731  xadccellx += cellADC*cellXYZ[0];
3732  zadccellz += cellADC*cellXYZ[2];
3733  }
3734 
3735  for(unsigned int nh=0; nh<cluster.NYCell();nh++){
3736  art::Ptr<rb::CellHit> cellhit = cluster.YCell(nh);
3737  Double_t cellXYZ[3]={0,0,0};
3739  geom->Plane(cellhit->Plane())->Cell(cellhit->Cell())->GetCenter(cellXYZ);
3740  double cellADC = double(cellhit->ADC());
3741  yadc += cellADC;
3742  zadc += cellADC;
3743  yadccelly += cellADC*cellXYZ[1];
3744  zadccellz += cellADC*cellXYZ[2];
3745  }
3746  if(xadc>0)cx = xadccellx/xadc;
3747  if(yadc>0)cy = yadccelly/yadc;
3748  if(zadc>0)cz = zadccellz/zadc;
3749  TVector3 centroid(cx, cy, cz);
3750  return centroid;
3751  }
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
unsigned short Plane() const
Definition: CellHit.h:39
const PlaneGeo * Plane(unsigned int i) const
unsigned short Cell() const
Definition: CellHit.h:40
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
void geom(int which=0)
Definition: geom.C:163
double jmshower::RecoJMShower::GetDedxInvLongLL ( rb::Prong  shower,
const int  type 
)

Definition at line 4594 of file RecoJMShower_module.cc.

References make_cached_def::e0, make_cached_def::e1, Energy(), fShowerPClusterCol, GetInterPlaneDedxProb(), GetPlaneDedx(), GetShwStop(), MECModelEnuComparisons::i, rb::Cluster::ID(), test_ParserArtEvents::log, pBadChannel, NDAPDHVSetting::plane, elec2geo::pos, gen_flatrecord::size, and rb::Prong::Start().

Referenced by RecoShowers().

4594  {
4595  int i = shower.ID();
4596  TVector3 pos((shower.Start()[0]+RecoJMShower::GetShwStop(shower)[0])/2.0,(shower.Start()[1]+RecoJMShower::GetShwStop(shower)[1])/2.0,(shower.Start()[2]+RecoJMShower::GetShwStop(shower)[2])/2.0);
4597  int iReg = 0;
4598  if(pos[0]>=0&&pos[1]>=0)iReg = 0;
4599  if(pos[0]>=0&&pos[1]<0)iReg = 1;
4600  if(pos[0]<0&&pos[1]>=0)iReg = 2;
4601  if(pos[0]<=0&&pos[1]<=0)iReg = 3;
4602 
4603  double gev = RecoJMShower::Energy(shower);
4604  double elower[11] = {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75};
4605  double eupper[11] = {0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25};
4606  double epoint[11];
4607  for(int ig=0;ig<11;ig++){
4608  epoint[ig] = (elower[ig]+eupper[ig])/2.;
4609  }
4610  int igev0 = 0;
4611  int igev1 = 0;
4612  if(gev>epoint[10]){igev0=10;igev1=10;}
4613  else if(gev<epoint[0]){igev0=0;igev1=0;}
4614  else{
4615  for(int ig=0;ig<10;ig++){
4616  if(gev>epoint[ig]&&gev<epoint[ig+1]){igev0=ig;igev1=ig+1;}
4617  }
4618  }
4619  if(igev0>10)igev0=10;
4620  if(igev1>10)igev1=10;
4621  double e0 = gev-epoint[igev0];
4622  double e1 = epoint[igev1] - gev;
4623 
4624 
4625  double sumlgprob = 0;
4626 // unsigned int totplane = 90;
4627  unsigned int totplane = fShowerPClusterCol[i].size();
4628  unsigned int nonzeroplane = 0;
4629 
4630 /*
4631  unsigned int shnp = fShowerPClusterCol[i].size();
4632  if(int(jmshower::kElectron)==type)totplane=std::max(shnp,fHtExpPlaneE[iReg][igev1]);
4633  if(int(jmshower::kGamma)==type)totplane=std::max(shnp,fHtExpPlaneG[iReg][igev1]);
4634  if(int(jmshower::kMuon)==type)totplane=std::max(shnp,fHtExpPlaneMu[iReg][igev1]);
4635  if(int(jmshower::kPi0)==type)totplane=std::max(shnp,fHtExpPlanePi0[iReg][igev1]);
4636  if(int(jmshower::kHad)==type)totplane=std::max(shnp,fHtExpPlaneHad[iReg][igev1]);
4637  if(int(jmshower::kProton)==type)totplane=std::max(shnp,fHtExpPlaneP[iReg][igev1]);
4638  if(int(jmshower::kNeutron)==type)totplane=std::max(shnp,fHtExpPlaneN[iReg][igev1]);
4639  if(int(jmshower::kPion)==type)totplane=std::max(shnp,fHtExpPlanePi[iReg][igev1]);
4640  if(int(jmshower::kElectronQE)==type)totplane=std::max(shnp,fHtExpPlaneEqe[iReg][igev1]);
4641  if(int(jmshower::kElectronRES)==type)totplane=std::max(shnp,fHtExpPlaneEres[iReg][igev1]);
4642  if(int(jmshower::kElectronDIS)==type)totplane=std::max(shnp,fHtExpPlaneEdis[iReg][igev1]);
4643  if(int(jmshower::kElectronCOH)==type)totplane=std::max(shnp,fHtExpPlaneEcoh[iReg][igev1]);
4644 */
4645 
4646 
4647  for(unsigned int plane=0;plane<totplane;plane++){
4648 // TVector3 disttoedge = RecoJMShower::GetTrkPlaneDistToEdge(fShowerPClusterCol[i][plane].Plane(), shower);
4649 // double planerad = RecoJMShower::Radius(fShowerPClusterCol[i][plane], shower);
4650 // if(disttoedge.x()<2*planerad)continue;
4651 // if(disttoedge.y()<2*planerad)continue;
4652 // if(disttoedge.z()<2*planerad)continue;
4653 // double path = RecoJMShower::GetTrkHitPath(fShowerPClusterCol[i][plane].Plane(), 0, shower);
4654 // double ep = RecoJMShower::Energy(fShowerPClusterCol[i][plane],shower.ID());
4655 // double dedx = ep/path;
4656  double dedx = 0;
4657  if(totplane-plane<fShowerPClusterCol[i].size()){
4658  dedx = RecoJMShower::GetPlaneDedx(shower, totplane-plane);
4659  }
4660  double prob = 0;
4661 // unsigned int nmip = RecoJMShower::NMIPPlane(shower);
4662  prob = RecoJMShower::GetInterPlaneDedxProb(shower, plane, dedx, type, pos, iReg, igev0, igev1, e0, e1);
4663 // prob = RecoJMShower::GetPlaneDedxProb(shower, i, hiteweightmap, 2,plane, dedx, type);
4664  if(prob<0.0001)prob=0.0001;
4665  double lgprob = log(prob);
4666  if(pBadChannel==true){ //Only consider non-zero plane
4667  if(dedx>0.0000001){
4668  nonzeroplane++;
4669  sumlgprob+=lgprob;
4670  }
4671  }else{
4672  nonzeroplane++;
4673  sumlgprob+=lgprob;
4674  }
4675  }
4676  if(pBadChannel==true){
4677  if(nonzeroplane>0)sumlgprob = sumlgprob/(double)nonzeroplane;
4678  }else{
4679  if(totplane>0)sumlgprob = sumlgprob/(double)totplane;
4680  }
4681 // sumlgprob = sumlgprob/double(fShowerPClusterCol[i].size());
4682  if(sumlgprob>999)sumlgprob=999;
4683  if(sumlgprob<-999)sumlgprob=-999;
4684  return sumlgprob;
4685  }
TVector3 GetShwStop(rb::Prong shw)
double Energy(rb::Cluster cluster, int i)
double GetInterPlaneDedxProb(rb::Prong shower, unsigned int plane, double dedx, const int type, TVector3 pos, int iReg, int igev0, int igev1, double e0, double e1)
virtual TVector3 Start() const
Definition: Prong.h:73
std::vector< std::vector< jmshower::PCluster > > fShowerPClusterCol
double GetPlaneDedx(rb::Prong shower, unsigned int plane)
const int ID() const
Definition: Cluster.h:75
double jmshower::RecoJMShower::GetDedxLongLL ( rb::Prong  shower,
const int  type 
)

Definition at line 4462 of file RecoJMShower_module.cc.

References make_cached_def::e0, make_cached_def::e1, Energy(), fShowerPClusterCol, GetInterPlaneDedxProb(), GetPlaneDedx(), GetShwStop(), MECModelEnuComparisons::i, rb::Cluster::ID(), jmshower::kElectronSG, test_ParserArtEvents::log, pBadChannel, NDAPDHVSetting::plane, elec2geo::pos, gen_flatrecord::size, and rb::Prong::Start().

Referenced by RecoShowers().

4462  {
4463  int i = shower.ID();
4464  TVector3 pos((shower.Start()[0]+RecoJMShower::GetShwStop(shower)[0])/2.0,(shower.Start()[1]+RecoJMShower::GetShwStop(shower)[1])/2.0,(shower.Start()[2]+RecoJMShower::GetShwStop(shower)[2])/2.0);
4465  int iReg = 0;
4466  if(pos[0]>=0&&pos[1]>=0)iReg = 0;
4467  if(pos[0]>=0&&pos[1]<0)iReg = 1;
4468  if(pos[0]<0&&pos[1]>=0)iReg = 2;
4469  if(pos[0]<=0&&pos[1]<=0)iReg = 3;
4470 
4471  double gev = RecoJMShower::Energy(shower);
4472  double elower[11] = {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75};
4473  double eupper[11] = {0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25};
4474  double epoint[11];
4475  for(int ig=0;ig<11;ig++){
4476  epoint[ig] = (elower[ig]+eupper[ig])/2.;
4477  }
4478  int igev0 = 0;
4479  int igev1 = 0;
4480  if(gev>epoint[10]){igev0=10;igev1=10;}
4481  else if(gev<epoint[0]){igev0=0;igev1=0;}
4482  else{
4483  for(int ig=0;ig<10;ig++){
4484  if(gev>epoint[ig]&&gev<epoint[ig+1]){igev0=ig;igev1=ig+1;}
4485  }
4486  }
4487  if(igev0>10)igev0=10;
4488  if(igev1>10)igev1=10;
4489  double e0 = gev-epoint[igev0];
4490  double e1 = epoint[igev1] - gev;
4491 
4492  double elowerm[41]= {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75,
4493  5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75,
4494  10.25,10.75,11.25,11.75,12.25,12.75,13.25,13.75,14.25,14.75,
4495  15.25,15.75,16.25,16.75,17.25,17.75,18.25,18.75,19.25,19.75
4496  };
4497 
4498  double eupperm[41]= { 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75,
4499  5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75,
4500  10.25,10.75,11.25,11.75,12.25,12.75,13.25,13.75,14.25,14.75,
4501  15.25,15.75,16.25,16.75,17.25,17.75,18.25,18.75,19.25,19.75, 20.25
4502  };
4503  double epointm[41];
4504  for(int ig=0;ig<41;ig++){
4505  epointm[ig] = (elowerm[ig]+eupperm[ig])/2.;
4506  }
4507 
4508  int igevm0 = 0;
4509  int igevm1 = 0;
4510  if(gev>epointm[40]){igevm0=40;igevm1=40;}
4511  else if(gev<epointm[0]){igevm0=0;igevm1=0;}
4512  else{
4513  for(int ig=0;ig<40;ig++){
4514  if(gev>epointm[ig]&&gev<epointm[ig+1]){igevm0=ig;igevm1=ig+1;}
4515  }
4516  }
4517  if(igevm0>40)igevm0=40;
4518  if(igevm1>40)igevm1=40;
4519  double em0 = gev-epointm[igevm0];
4520  double em1 = epointm[igevm1] - gev;
4521 
4522 
4523  double sumlgprob = 0;
4524 // unsigned int totplane = 90;
4525  unsigned int totplane = fShowerPClusterCol[i].size();
4526  unsigned int nonzeroplane = 0;
4527 
4528 /*
4529  unsigned int shnp = fShowerPClusterCol[i].size();
4530  if(int(jmshower::kElectron)==type)totplane=std::max(shnp,fHtExpPlaneE[iReg][igev1]);
4531  if(int(jmshower::kGamma)==type)totplane=std::max(shnp,fHtExpPlaneG[iReg][igev1]);
4532  if(int(jmshower::kMuon)==type)totplane=std::max(shnp,fHtExpPlaneMu[iReg][igev1]);
4533  if(int(jmshower::kPi0)==type)totplane=std::max(shnp,fHtExpPlanePi0[iReg][igev1]);
4534  if(int(jmshower::kHad)==type)totplane=std::max(shnp,fHtExpPlaneHad[iReg][igev1]);
4535  if(int(jmshower::kProton)==type)totplane=std::max(shnp,fHtExpPlaneP[iReg][igev1]);
4536  if(int(jmshower::kNeutron)==type)totplane=std::max(shnp,fHtExpPlaneN[iReg][igev1]);
4537  if(int(jmshower::kPion)==type)totplane=std::max(shnp,fHtExpPlanePi[iReg][igev1]);
4538  if(int(jmshower::kElectronQE)==type)totplane=std::max(shnp,fHtExpPlaneEqe[iReg][igev1]);
4539  if(int(jmshower::kElectronRES)==type)totplane=std::max(shnp,fHtExpPlaneEres[iReg][igev1]);
4540  if(int(jmshower::kElectronDIS)==type)totplane=std::max(shnp,fHtExpPlaneEdis[iReg][igev1]);
4541  if(int(jmshower::kElectronCOH)==type)totplane=std::max(shnp,fHtExpPlaneEcoh[iReg][igev1]);
4542 */
4543 
4544 
4545  //std::cout<<"type, energy, cos, np, totplane "<<type<<" "<<gev<<" "<<shower.Dir()[2]<<" "<<shnp<<" "<<totplane<<std::endl;
4546 
4547  for(unsigned int plane=0;plane<totplane;plane++){
4548 // TVector3 disttoedge = RecoJMShower::GetTrkPlaneDistToEdge(fShowerPClusterCol[i][plane].Plane(), shower);
4549 // double planerad = RecoJMShower::Radius(fShowerPClusterCol[i][plane], shower);
4550 // if(disttoedge.x()<2*planerad)continue;
4551 // if(disttoedge.y()<2*planerad)continue;
4552 // if(disttoedge.z()<2*planerad)continue;
4553 // double path = RecoJMShower::GetTrkHitPath(fShowerPClusterCol[i][plane].Plane(), 0, shower);
4554 // double ep = RecoJMShower::Energy(fShowerPClusterCol[i][plane],shower.ID());
4555 // double dedx = ep/path;
4556  double dedx = 0;
4557  if(plane<fShowerPClusterCol[i].size()){
4558  dedx = RecoJMShower::GetPlaneDedx(shower, plane);
4559  }
4560  double prob = 0;
4561 // unsigned int nmip = RecoJMShower::NMIPPlane(shower);
4562  if(type != int(jmshower::kElectronSG)){
4563  prob = RecoJMShower::GetInterPlaneDedxProb(shower, plane, dedx, type, pos, iReg, igev0, igev1, e0, e1);
4564  }else{
4565  prob = RecoJMShower::GetInterPlaneDedxProb(shower, plane, dedx, type, pos, iReg, igevm0, igevm1, em0, em1);
4566  }
4567 // prob = RecoJMShower::GetPlaneDedxProb(shower, i, hiteweightmap, 2,plane, dedx, type);
4568  if(prob<0.0001)prob=0.0001;
4569  double lgprob = log(prob);
4570  if(pBadChannel==true){ //Only consider non-zero plane
4571  if(dedx>0.0000001){
4572  nonzeroplane++;
4573  sumlgprob+=lgprob;
4574  }
4575  }else{
4576  nonzeroplane++;
4577  sumlgprob+=lgprob;
4578  }
4579  //std::cout<<"type: "<<type<<" "<<plane<<" "<<prob<<" "<<sumlgprob<<std::endl;
4580  }
4581  if(pBadChannel==true){
4582  if(nonzeroplane>0)sumlgprob = sumlgprob/(double)nonzeroplane;
4583  }else{
4584  if(totplane>0)sumlgprob = sumlgprob/(double)totplane;
4585  }
4586 // std::cout<<"sumlgprob, totplane "<<sumlgprob<<" "<<totplane<<std::endl;
4587 // sumlgprob = sumlgprob/double(fShowerPClusterCol[i].size());
4588  if(sumlgprob>999)sumlgprob=999;
4589  if(sumlgprob<-999)sumlgprob=-999;
4590  //std::cout<<"DedxLongLL type: "<<type<<" ll "<<sumlgprob<<std::endl;
4591  return sumlgprob;
4592  }
TVector3 GetShwStop(rb::Prong shw)
double Energy(rb::Cluster cluster, int i)
double GetInterPlaneDedxProb(rb::Prong shower, unsigned int plane, double dedx, const int type, TVector3 pos, int iReg, int igev0, int igev1, double e0, double e1)
virtual TVector3 Start() const
Definition: Prong.h:73
std::vector< std::vector< jmshower::PCluster > > fShowerPClusterCol
double GetPlaneDedx(rb::Prong shower, unsigned int plane)
const int ID() const
Definition: Cluster.h:75
double jmshower::RecoJMShower::GetDedxLongLL ( rb::Prong  shower,
const int  type,
unsigned int  startp,
int  contp 
)
double jmshower::RecoJMShower::GetDedxTransLL ( rb::Prong  shower,
const int  type 
)

Definition at line 5087 of file RecoJMShower_module.cc.

References make_cached_def::e0, make_cached_def::e1, Energy(), GetCellTransDedx(), GetInterCellTransDedxProb(), GetShwStop(), jmshower::kElectronSG, test_ParserArtEvents::log, elec2geo::pos, and rb::Prong::Start().

Referenced by RecoShowers().

5087  {
5088 // int i = shower.ID();
5089  double transcelldedx[20];
5090  for(int itd=0;itd<20;itd++){
5091  transcelldedx[itd] = 0;
5092  }
5093  for(unsigned itp = 0;itp<20;itp++){
5094  transcelldedx[itp] = RecoJMShower::GetCellTransDedx(shower, itp);
5095  }
5096  double sumlgprob = 0;
5097 
5098  TVector3 pos((shower.Start()[0]+RecoJMShower::GetShwStop(shower)[0])/2.0,(shower.Start()[1]+RecoJMShower::GetShwStop(shower)[1])/2.0,(shower.Start()[2]+RecoJMShower::GetShwStop(shower)[2])/2.0);
5099  int iReg = 0;
5100  if(pos[0]>=0&&pos[1]>=0)iReg = 0;
5101  if(pos[0]>=0&&pos[1]<0)iReg = 1;
5102  if(pos[0]<0&&pos[1]>=0)iReg = 2;
5103  if(pos[0]<=0&&pos[1]<=0)iReg = 3;
5104 
5105 
5106  double gev = RecoJMShower::Energy(shower);
5107  double elower[11] = {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75};
5108  double eupper[11] = {0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25};
5109  double epoint[11];
5110  for(int ig=0;ig<11;ig++){
5111  epoint[ig] = (elower[ig]+eupper[ig])/2.;
5112  }
5113  int igev0 = 0;
5114  int igev1 = 0;
5115  if(gev>epoint[10]){igev0=10;igev1=10;}
5116  else if(gev<epoint[0]){igev0=0;igev1=0;}
5117  else{
5118  for(int ig=0;ig<10;ig++){
5119  if(gev>epoint[ig]&&gev<epoint[ig+1]){igev0=ig;igev1=ig+1;}
5120  }
5121  }
5122  if(igev0>10)igev0=10;
5123  if(igev1>10)igev1=10;
5124 
5125  double e0 = gev-epoint[igev0];
5126  double e1 = epoint[igev1] - gev;
5127 
5128 
5129  double elowerm[41]= {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75,
5130  5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75,
5131  10.25,10.75,11.25,11.75,12.25,12.75,13.25,13.75,14.25,14.75,
5132  15.25,15.75,16.25,16.75,17.25,17.75,18.25,18.75,19.25,19.75
5133  };
5134 
5135  double eupperm[41]= { 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75,
5136  5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75,
5137  10.25,10.75,11.25,11.75,12.25,12.75,13.25,13.75,14.25,14.75,
5138  15.25,15.75,16.25,16.75,17.25,17.75,18.25,18.75,19.25,19.75, 20.25
5139  };
5140  double epointm[41];
5141  for(int ig=0;ig<41;ig++){
5142  epointm[ig] = (elowerm[ig]+eupperm[ig])/2.;
5143  }
5144 
5145  int igevm0 = 0;
5146  int igevm1 = 0;
5147  if(gev>epointm[40]){igevm0=40;igevm1=40;}
5148  else if(gev<epointm[0]){igevm0=0;igevm1=0;}
5149  else{
5150  for(int ig=0;ig<40;ig++){
5151  if(gev>epointm[ig]&&gev<epointm[ig+1]){igevm0=ig;igevm1=ig+1;}
5152  }
5153  }
5154  if(igevm0>40)igevm0=40;
5155  if(igevm1>40)igevm1=40;
5156 
5157  double em0 = gev-epointm[igevm0];
5158  double em1 = epointm[igevm1] - gev;
5159 
5160  int ntcell = 0;
5161  for(int itd=0;itd<20;itd++){
5162  if(transcelldedx[itd]<-0.5)continue;
5163  double prob = 0;
5164  if(type != int(jmshower::kElectronSG)){
5165  prob = RecoJMShower::GetInterCellTransDedxProb(shower, itd, transcelldedx[itd], type, pos, iReg, igev0, igev1, e0, e1);
5166  }else{
5167  prob = RecoJMShower::GetInterCellTransDedxProb(shower, itd, transcelldedx[itd], type, pos, iReg, igevm0, igevm1, em0, em1);
5168  }
5169  if(prob<0.0001)prob=0.0001;
5170  double lgprob = log(prob);
5171  sumlgprob+=lgprob;
5172  ntcell++;
5173  }
5174  if(ntcell>0)sumlgprob = sumlgprob/(double)ntcell;
5175 // if(ntcell>0)sumlgprob = sumlgprob/20.;
5176  if(sumlgprob>999)sumlgprob=999;
5177  if(sumlgprob<-999)sumlgprob=-999;
5178  //std::cout<<"trans ll prob "<<sumlgprob<<std::endl;
5179  return sumlgprob;
5180  }
TVector3 GetShwStop(rb::Prong shw)
double Energy(rb::Cluster cluster, int i)
virtual TVector3 Start() const
Definition: Prong.h:73
double GetInterCellTransDedxProb(rb::Prong shower, unsigned int cell, double dedx, const int type, TVector3 pos, int iReg, int igev0, int igev1, double e0, double e1)
double GetCellTransDedx(rb::Prong shower, unsigned int cell)
double jmshower::RecoJMShower::GetInterCellTransDedxProb ( rb::Prong  shower,
unsigned int  cell,
double  dedx,
const int  type,
TVector3  pos,
int  iReg,
int  igev0,
int  igev1,
double  e0,
double  e1 
)

Definition at line 5037 of file RecoJMShower_module.cc.

References geo::GeometryBase::DetHalfHeight(), geo::GeometryBase::DetHalfWidth(), geom(), GetCellTransDedxProb(), X1, and Y1.

Referenced by GetDedxTransLL().

5037  {
5038 /*
5039  TVector3 pos((shower.Start()[0]+RecoJMShower::GetShwStop(shower)[0])/2.0,(shower.Start()[1]+RecoJMShower::GetShwStop(shower)[1])/2.0,(shower.Start()[2]+RecoJMShower::GetShwStop(shower)[2])/2.0);
5040  int ireg = 0;
5041  if(pos[0]>=0&&pos[1]>=0)ireg = 0;
5042  if(pos[0]>=0&&pos[1]<0)ireg = 1;
5043  if(pos[0]<0&&pos[1]>=0)ireg = 2;
5044  if(pos[0]<=0&&pos[1]<=0)ireg = 3;
5045 */
5047 
5048  double probPXPY = 0;
5049  bool interPos = false;
5050 
5051 
5052  if(interPos==true){
5053  double X0 = -geom->DetHalfWidth()/2.0;
5054  double X1 = geom->DetHalfWidth()/2.0;
5055  double Y0 = -geom->DetHalfHeight()/2.0;
5056  double Y1 = geom->DetHalfHeight()/2.0;
5057 
5058  double probreg[4]={1,1,1,1};
5059  for(int ireg=0;ireg<4;ireg++){
5060  double probe0 = RecoJMShower::GetCellTransDedxProb(shower, ireg, igev0, cell, dedx, type);
5061  double probe1 = RecoJMShower::GetCellTransDedxProb(shower, ireg, igev1, cell, dedx, type);
5062  if(igev0!=igev1)probreg[ireg] = (e1*probe0+e0*probe1)/(e0+e1);
5063  if(igev0==igev1)probreg[ireg] = probe0;
5064  }
5065 
5066  double PX = pos[0];
5067  double PY = pos[1];
5068  double probX1Y1 = probreg[0];
5069  double probX1Y0 = probreg[1];
5070  double probX0Y1 = probreg[2];
5071  double probX0Y0 = probreg[3];
5072  double probX1PY = ((PY-Y0)/(Y1-Y0))*(probX1Y1-probX1Y0)+probX1Y0;
5073  double probX0PY = ((PY-Y0)/(Y1-Y0))*(probX0Y1-probX0Y0)+probX0Y0;
5074  probPXPY = ((PX-X0)/(X1-X0))*(probX1PY-probX0PY)+probX0PY;
5075  }
5076 
5077  if(interPos==false||probPXPY<0){
5078  double probe0 = RecoJMShower::GetCellTransDedxProb(shower, iReg, igev0, cell, dedx, type);
5079  double probe1 = RecoJMShower::GetCellTransDedxProb(shower, iReg, igev1, cell, dedx, type);
5080  if(igev0!=igev1)probPXPY = (e1*probe0+e0*probe1)/(e0+e1);
5081  if(igev0==igev1)probPXPY = probe0;
5082  }
5083 
5084  return probPXPY;
5085  }
Double_t X1
Definition: plot.C:264
Double_t Y1
Definition: plot.C:264
double DetHalfHeight() const
double DetHalfWidth() const
void geom(int which=0)
Definition: geom.C:163
double GetCellTransDedxProb(rb::Prong shower, int ireg, int igev, unsigned int cell, double dedx, const int type)
double jmshower::RecoJMShower::GetInterPlaneDedxProb ( rb::Prong  shower,
unsigned int  plane,
double  dedx,
const int  type,
TVector3  pos,
int  iReg,
int  igev0,
int  igev1,
double  e0,
double  e1 
)

Definition at line 4310 of file RecoJMShower_module.cc.

References std::abs(), std::cos(), geo::GeometryBase::DetHalfHeight(), geo::GeometryBase::DetHalfWidth(), rb::Prong::Dir(), geom(), GetPlaneDedxProb(), makeTrainCVSamples::int, X1, and Y1.

Referenced by GetDedxInvLongLL(), GetDedxLongLL(), and RecoShowers().

4310  {
4311 
4312  double probPXPY = 0;
4313  bool interPos = false;
4314 
4315 // int i = shower.ID();
4317 // TVector3 pos = RecoJMShower::GetTrkPlanePos(fShowerPClusterCol[i][plane].Plane(), shower);
4318 // std::cout<<"in getinterplane pos X,Y "<<pos[0]<<" "<<pos[1]<<std::endl;
4319 
4320 // TVector3 pos((shower.Start()[0]+RecoJMShower::GetShwStop(shower)[0])/2.0,(shower.Start()[1]+RecoJMShower::GetShwStop(shower)[1])/2.0,(shower.Start()[2]+RecoJMShower::GetShwStop(shower)[2])/2.0);
4321 
4322 
4323  double cos = shower.Dir()[2];
4324 
4325  if(interPos==true){
4326  double X0 = -geom->DetHalfWidth()/2.0;
4327  double X1 = geom->DetHalfWidth()/2.0;
4328  double Y0 = -geom->DetHalfHeight()/2.0;
4329  double Y1 = geom->DetHalfHeight()/2.0;
4330 
4331  double probreg[4]={1,1,1,1};
4332  for(int ireg=0;ireg<4;ireg++){
4333  unsigned int plane0 = int(int(plane)/std::abs(cos));
4334  unsigned int plane1 = int(int(plane)/std::abs(cos)) + 1;
4335  double r0 = double(plane/std::abs(cos))-double(plane0);
4336  double r1 = 1-r0;
4337  double probe0r0 = 1;
4338  double probe0r1 = 1;
4339  double probe1r0 = 1;
4340  double probe1r1 = 1;
4341 
4342  probe0r0 = RecoJMShower::GetPlaneDedxProb(shower, ireg, igev0, plane0, dedx, type);
4343  probe0r1 = RecoJMShower::GetPlaneDedxProb(shower, ireg, igev0, plane1, dedx, type);
4344  probe1r0 = RecoJMShower::GetPlaneDedxProb(shower, ireg, igev1, plane0, dedx, type);
4345  probe1r1 = RecoJMShower::GetPlaneDedxProb(shower, ireg, igev1, plane1, dedx, type);
4346 // std::cout<<"ireg, prob probe0r0 "<<ireg<<" "<<probe0r0<<std::endl;
4347 // std::cout<<"ireg, prob probe0r1 "<<ireg<<" "<<probe0r1<<std::endl;
4348 // std::cout<<"ireg, prob probe1r0 "<<ireg<<" "<<probe1r0<<std::endl;
4349 // std::cout<<"ireg, prob probe1r0 "<<ireg<<" "<<probe1r0<<std::endl;
4350  double probe0 = (r1*probe0r0+r0*probe0r1);
4351  double probe1 = (r1*probe1r0+r0*probe1r1);
4352  if(igev0!=igev1)probreg[ireg] = (e1*probe0+e0*probe1)/(e0+e1);
4353  if(igev0==igev1)probreg[ireg] = probe0;
4354 // std::cout<<"ireg, probreg "<<ireg<<" "<<probreg[ireg]<<std::endl;
4355  }
4356  double PX = pos[0];
4357  double PY = pos[1];
4358  double probX1Y1 = probreg[0];
4359  double probX1Y0 = probreg[1];
4360  double probX0Y1 = probreg[2];
4361  double probX0Y0 = probreg[3];
4362 // std::cout<<"prob X1Y1, X1Y0, X0Y1, X0Y0 "<<probX1Y1<<" "<<probX1Y0<<" "<<probX0Y1<<" "<<probX0Y0<<std::endl;
4363  double probX1PY = ((PY-Y0)/(Y1-Y0))*(probX1Y1-probX1Y0)+probX1Y0;
4364  double probX0PY = ((PY-Y0)/(Y1-Y0))*(probX0Y1-probX0Y0)+probX0Y0;
4365  probPXPY = ((PX-X0)/(X1-X0))*(probX1PY-probX0PY)+probX0PY;
4366  }
4367 // std::cout<<"prob X1PY, X0PY, PXPY"<<probX1PY<<" "<<probX0PY<<" "<<probPXPY<<" "<<std::endl;
4368 // std::cout<<"gev,e0,e1,r0,r0: "<<gev<<" "<<e0<<" "<<e1<<" "<<r0<<" "<<r1<<std::endl;
4369 // std::cout<<"prob,pe0r0,pe0r1,pe1r0,pe1r1: "<<prob<<" "<<probe0r0<<" "<<probe0r1<<" "<<probe1r0<<" "<<probe1r1<<std::endl;
4370  if(interPos==false||probPXPY<0){
4371  unsigned int plane0 = int(int(plane)/std::abs(cos));
4372  unsigned int plane1 = int(int(plane)/std::abs(cos)) + 1;
4373  double r0 = double(plane/std::abs(cos))-double(plane0);
4374  double r1 = 1-r0;
4375  double probe0r0 = 1;
4376  double probe0r1 = 1;
4377  double probe1r0 = 1;
4378  double probe1r1 = 1;
4379  probe0r0 = RecoJMShower::GetPlaneDedxProb(shower, iReg, igev0, plane0, dedx, type);
4380  probe0r1 = RecoJMShower::GetPlaneDedxProb(shower, iReg, igev0, plane1, dedx, type);
4381  probe1r0 = RecoJMShower::GetPlaneDedxProb(shower, iReg, igev1, plane0, dedx, type);
4382  probe1r1 = RecoJMShower::GetPlaneDedxProb(shower, iReg, igev1, plane1, dedx, type);
4383  double probe0 = (r1*probe0r0+r0*probe0r1);
4384  double probe1 = (r1*probe1r0+r0*probe1r1);
4385  if(igev0!=igev1)probPXPY = (e1*probe0+e0*probe1)/(e0+e1);
4386  if(igev0==igev1)probPXPY = probe0;
4387  }
4388 
4389  return probPXPY;
4390  }
double GetPlaneDedxProb(rb::Prong shower, int ireg, int igev, unsigned int plane, double dedx, const int type)
float abs(float number)
Definition: d0nt_math.hpp:39
Double_t X1
Definition: plot.C:264
Double_t Y1
Definition: plot.C:264
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
double DetHalfHeight() const
double DetHalfWidth() const
void geom(int which=0)
Definition: geom.C:163
T cos(T number)
Definition: d0nt_math.hpp:78
int jmshower::RecoJMShower::GetPlaneCentroidCell ( jmshower::PCluster  pcluster,
rb::Prong  shower 
)

Definition at line 5811 of file RecoJMShower_module.cc.

References rb::CellHit::Cell(), rb::Cluster::Cell(), DEFINE_ART_MODULE(), makeTrainCVSamples::int, rb::Cluster::NCell(), and rb::CellHit::PE().

Referenced by RecoShowers().

5811  {
5812  double gev = 0;
5813  double gevcell = 0;
5814  if(pcluster.NCell()==0)return -1;
5815  for(unsigned int nh=0; nh<pcluster.NCell();nh++){
5816  art::Ptr<rb::CellHit> cellhit = pcluster.Cell(nh);
5817  double cellgev = cellhit->PE();
5818  gev += cellgev;
5819  gevcell += cellgev*(double)cellhit->Cell();
5820  }
5821  int planecentroidcell = (int)pcluster.Cell(0)->Cell();
5822  if(gev>0.000001) planecentroidcell = (int)gevcell/gev;
5823  return planecentroidcell;
5824  }
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
unsigned short Cell() const
Definition: CellHit.h:40
float PE() const
Definition: CellHit.h:42
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 jmshower::RecoJMShower::GetPlaneDedx ( rb::Prong  shower,
unsigned int  plane 
)

Definition at line 3760 of file RecoJMShower_module.cc.

References om::cout, allTimeWatchdog::endl, Energy(), lem_server::ep, fDetid, fShowerPClusterCol, GetTrkHitPath(), MECModelEnuComparisons::i, rb::Cluster::ID(), isrealdata, path, pCorDataDedx, PlaneDedxDataMC(), and gen_flatrecord::size.

Referenced by CalcANN(), ContStartPlane(), GetDedxInvLongLL(), GetDedxLongLL(), NMIPPlane(), and RecoShowers().

3760  {
3761  int i = shower.ID();
3762  if(fShowerPClusterCol[i].size()<1)
3763  {
3764  std::cout<<"JMShower: Nothing in plane "<<i<<std::endl;
3765  return 9999;
3766  }
3767  double path = RecoJMShower::GetTrkHitPath(fShowerPClusterCol[i][plane].Plane(), 0, shower);
3768  double ep = RecoJMShower::Energy(fShowerPClusterCol[i][plane], shower.ID());
3769 // TVector3 hitpos = RecoJMShower::GetTrkPlanePos(fShowerPClusterCol[i][plane].Plane(), shower);
3770 // double poscor = (0.927377+6.46394e-05*hitpos[0])*(0.993889+3.83504e-05*hitpos[1]);/from plane
3771 // double poscor = (9.9899+7.46983e-05*avgx+7.34449e-05*avgy);//from shower
3772 // double dedx = (ep/poscor)/path;
3773  double dedx = ep/path;
3774  if(isrealdata==true&&pCorDataDedx)dedx=dedx/RecoJMShower::PlaneDedxDataMC(dedx, fDetid);
3775  //std::cout<<"shower: "<<i<<" plane idx: "<<plane<<" plane: "<<fShowerPClusterCol[i][plane].Plane()<<" path: "<<path<<" ep: "<<ep<<" dedx: "<<dedx<<std::endl;
3776  return dedx;
3777  }
double Energy(rb::Cluster cluster, int i)
std::vector< std::vector< jmshower::PCluster > > fShowerPClusterCol
double PlaneDedxDataMC(double gev, int detid)
double GetTrkHitPath(unsigned int plane, unsigned int cell, rb::Prong shower)
OStream cout
Definition: OStream.cxx:6
const std::string path
Definition: plot_BEN.C:43
const int ID() const
Definition: Cluster.h:75
double jmshower::RecoJMShower::GetPlaneDedxProb ( rb::Prong  shower,
int  ireg,
int  igev,
unsigned int  plane,
double  dedx,
const int  type 
)

Definition at line 3780 of file RecoJMShower_module.cc.

References fHtPlaneDedxE, fHtPlaneDedxEcoh, fHtPlaneDedxEdis, fHtPlaneDedxEqe, fHtPlaneDedxEres, fHtPlaneDedxEsg, fHtPlaneDedxG, fHtPlaneDedxMu, fHtPlaneDedxN, fHtPlaneDedxP, fHtPlaneDedxPi, fHtPlaneDedxPi0, GetEntries(), HandyFuncs::Interpolate(), jmshower::kElectron, jmshower::kElectronCOH, jmshower::kElectronDIS, jmshower::kElectronQE, jmshower::kElectronRES, jmshower::kElectronSG, jmshower::kGamma, jmshower::kMuon, jmshower::kNeutron, jmshower::kPi0, jmshower::kPion, jmshower::kProton, getGoodRuns4SAM::n, and NDAPDHVSetting::plane.

Referenced by GetInterPlaneDedxProb().

3780  {
3781 // int i = shower.ID();
3782  double p = 0;
3783 
3784  if(int(jmshower::kElectronSG)==type){
3785  if(ireg<0||ireg>3)return 1;
3786  if(igev>=41||igev<0)return 1;
3787  if(plane>=200)return 1;
3788  }else{
3789  if(ireg<0||ireg>3)return 1;
3790  if(igev>=11||igev<0)return 1;
3791  if(plane>=200)return 1;
3792  }
3793 
3794  if(int(jmshower::kElectron)==type){
3795  double ntot = double(fHtPlaneDedxE[ireg][igev][plane]->GetEntries());
3796  double n = double(fHtPlaneDedxE[ireg][igev][plane]->Interpolate(dedx));
3797  double nbin = fHtPlaneDedxE[ireg][igev][plane]->GetNbinsX();
3798  if(n<0.0001)n=0.0001;
3799  if(ntot<1){
3800  p = 0;
3801  }else{
3802  p = n*nbin/ntot;
3803  }
3804  //std::cout<<"type "<<type<<" "<<ireg<<" "<<igev<<" "<<plane<<" "<<dedx<<" "<<ntot<<" "<<n<<" "<<nbin<<" "<<p<<std::endl;
3805  }
3806  if(int(jmshower::kGamma)==type){
3807  double ntot = double(fHtPlaneDedxG[ireg][igev][plane]->GetEntries());
3808  double n = double(fHtPlaneDedxG[ireg][igev][plane]->Interpolate(dedx));
3809  double nbin = fHtPlaneDedxG[ireg][igev][plane]->GetNbinsX();
3810  if(n<0.0001)n=0.0001;
3811  if(ntot<1){
3812  p = 0;
3813  }else{
3814  p = n*nbin/ntot;
3815  }
3816  }
3817  if(int(jmshower::kMuon)==type){
3818  double ntot = double(fHtPlaneDedxMu[ireg][igev][plane]->GetEntries());
3819  double n = double(fHtPlaneDedxMu[ireg][igev][plane]->Interpolate(dedx));
3820  double nbin = fHtPlaneDedxMu[ireg][igev][plane]->GetNbinsX();
3821  if(n<0.0001)n=0.0001;
3822  if(ntot<1){
3823  p = 0;
3824  }else{
3825  p = n*nbin/ntot;
3826  }
3827  }
3828  if(int(jmshower::kPi0)==type){
3829  double ntot = double(fHtPlaneDedxPi0[ireg][igev][plane]->GetEntries());
3830  double n = double(fHtPlaneDedxPi0[ireg][igev][plane]->Interpolate(dedx));
3831  double nbin = fHtPlaneDedxPi0[ireg][igev][plane]->GetNbinsX();
3832  if(n<0.0001)n=0.0001;
3833  if(ntot<1){
3834  p = 0;
3835  }else{
3836  p = n*nbin/ntot;
3837  }
3838  }
3839 /*
3840  if(int(jmshower::kHad)==type){
3841  double ntot = double(fHtPlaneDedxHad[ireg][igev][plane]->GetEntries());
3842  double n = double(fHtPlaneDedxHad[ireg][igev][plane]->Interpolate(dedx));
3843  double nbin = fHtPlaneDedxHad[ireg][igev][plane]->GetNbinsX();
3844  if(n<0.0001)n=0.0001;
3845  if(ntot<1){
3846  p = 0;
3847  }else{
3848  p = n*nbin/ntot;
3849  }
3850  }
3851 */
3852 
3853  if(int(jmshower::kProton)==type){
3854  double ntot = double(fHtPlaneDedxP[ireg][igev][plane]->GetEntries());
3855  double n = double(fHtPlaneDedxP[ireg][igev][plane]->Interpolate(dedx));
3856  double nbin = fHtPlaneDedxP[ireg][igev][plane]->GetNbinsX();
3857  if(n<0.0001)n=0.0001;
3858  if(ntot<1){
3859  p = 0;
3860  }else{
3861  p = n*nbin/ntot;
3862  }
3863  }
3864 
3865  if(int(jmshower::kNeutron)==type){
3866  double ntot = double(fHtPlaneDedxN[ireg][igev][plane]->GetEntries());
3867  double n = double(fHtPlaneDedxN[ireg][igev][plane]->Interpolate(dedx));
3868  double nbin = fHtPlaneDedxN[ireg][igev][plane]->GetNbinsX();
3869  if(n<0.0001)n=0.0001;
3870  if(ntot<1){
3871  p = 0;
3872  }else{
3873  p = n*nbin/ntot;
3874  }
3875  }
3876 
3877  if(int(jmshower::kPion)==type){
3878  double ntot = double(fHtPlaneDedxPi[ireg][igev][plane]->GetEntries());
3879  double n = double(fHtPlaneDedxPi[ireg][igev][plane]->Interpolate(dedx));
3880  double nbin = fHtPlaneDedxPi[ireg][igev][plane]->GetNbinsX();
3881  if(n<0.0001)n=0.0001;
3882  if(ntot<1){
3883  p = 0;
3884  }else{
3885  p = n*nbin/ntot;
3886  }
3887  }
3888  if(int(jmshower::kElectronQE)==type){
3889  double ntot = double(fHtPlaneDedxEqe[ireg][igev][plane]->GetEntries());
3890  double n = double(fHtPlaneDedxEqe[ireg][igev][plane]->Interpolate(dedx));
3891  double nbin = fHtPlaneDedxEqe[ireg][igev][plane]->GetNbinsX();
3892  if(n<0.0001)n=0.0001;
3893  if(ntot<1){
3894  p = 0;
3895  }else{
3896  p = n*nbin/ntot;
3897  }
3898  }
3899 
3900  if(int(jmshower::kElectronRES)==type){
3901  double ntot = double(fHtPlaneDedxEres[ireg][igev][plane]->GetEntries());
3902  double n = double(fHtPlaneDedxEres[ireg][igev][plane]->Interpolate(dedx));
3903  double nbin = fHtPlaneDedxEres[ireg][igev][plane]->GetNbinsX();
3904  if(n<0.0001)n=0.0001;
3905  if(ntot<1){
3906  p = 0;
3907  }else{
3908  p = n*nbin/ntot;
3909  }
3910  }
3911 
3912  if(int(jmshower::kElectronDIS)==type){
3913  double ntot = double(fHtPlaneDedxEdis[ireg][igev][plane]->GetEntries());
3914  double n = double(fHtPlaneDedxEdis[ireg][igev][plane]->Interpolate(dedx));
3915  double nbin = fHtPlaneDedxEdis[ireg][igev][plane]->GetNbinsX();
3916  if(n<0.0001)n=0.0001;
3917  if(ntot<1){
3918  p = 0;
3919  }else{
3920  p = n*nbin/ntot;
3921  }
3922  }
3923 
3924  if(int(jmshower::kElectronCOH)==type){
3925  double ntot = double(fHtPlaneDedxEcoh[ireg][igev][plane]->GetEntries());
3926  double n = double(fHtPlaneDedxEcoh[ireg][igev][plane]->Interpolate(dedx));
3927  double nbin = fHtPlaneDedxEcoh[ireg][igev][plane]->GetNbinsX();
3928  if(n<0.0001)n=0.0001;
3929  if(ntot<1){
3930  p = 0;
3931  }else{
3932  p = n*nbin/ntot;
3933  }
3934  }
3935 
3936  if(int(jmshower::kElectronSG)==type){
3937  double ntot = double(fHtPlaneDedxEsg[ireg][igev][plane]->GetEntries());
3938  double n = double(fHtPlaneDedxEsg[ireg][igev][plane]->Interpolate(dedx));
3939  double nbin = fHtPlaneDedxEsg[ireg][igev][plane]->GetNbinsX();
3940  if(n<0.0001)n=0.0001;
3941  if(ntot<1){
3942  p = 0;
3943  }else{
3944  p = n*nbin/ntot;
3945  }
3946  }
3947  return p;
3948  }
TH1D * fHtPlaneDedxEres[4][11][200]
TH1D * fHtPlaneDedxPi0[4][11][200]
TH1D * fHtPlaneDedxMu[4][11][200]
TH1D * fHtPlaneDedxN[4][11][200]
const char * p
Definition: xmltok.h:285
TH1D * fHtPlaneDedxE[4][11][200]
TH1D * fHtPlaneDedxG[4][11][200]
cout<< t1-> GetEntries()
Definition: plottest35.C:29
TH1D * fHtPlaneDedxEcoh[4][11][200]
TH1D * fHtPlaneDedxEsg[4][41][200]
TH1D * fHtPlaneDedxEdis[4][11][200]
def Interpolate(x1, y1, x2, y2, yvalue)
Definition: HandyFuncs.py:218
TH1D * fHtPlaneDedxEqe[4][11][200]
TH1D * fHtPlaneDedxP[4][11][200]
TH1D * fHtPlaneDedxPi[4][11][200]
int jmshower::RecoJMShower::GetPlaneE1Cell ( jmshower::PCluster  pcluster)

Definition at line 5796 of file RecoJMShower_module.cc.

References rb::CellHit::Cell(), rb::Cluster::Cell(), makeTrainCVSamples::int, rb::Cluster::NCell(), and rb::CellHit::PE().

Referenced by RecoShowers().

5796  {
5797  double e1gev = 0;
5798  int e1gevcell = -1;
5799  if(pcluster.NCell()==0)return -1;
5800  for(unsigned int nh=0; nh<pcluster.NCell();nh++){
5801  art::Ptr<rb::CellHit> cellhit = pcluster.Cell(nh);
5802  double cellgev = cellhit->PE();
5803  if(cellgev>e1gev){
5804  e1gev = cellgev;
5805  e1gevcell = (int)cellhit->Cell();
5806  }
5807  }
5808  return e1gevcell;
5809  }
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
unsigned short Cell() const
Definition: CellHit.h:40
float PE() const
Definition: CellHit.h:42
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 jmshower::RecoJMShower::GetPointDoca ( TVector3  vtx,
TVector3  start,
TVector3  stop 
)

Definition at line 2908 of file RecoJMShower_module.cc.

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

Referenced by RecoShowers().

2908  {
2909  double d = 9999;
2910  double denorm = (start-stop).Mag();
2911  if(denorm<0.00001){d = (vtx-start).Mag();}else{
2912  d = ((vtx-start).Cross(vtx-stop)).Mag()/denorm;
2913  TVector3 dir = stop-start;
2914  }
2915  return d;
2916  }
Float_t d
Definition: plot.C:236
TDirectory * dir
Definition: macro.C:5
float Mag() const
template<typename PROD , BranchType B>
ProductID art::EDProducer::getProductID ( std::string const &  instanceName = {}) const
inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

Referenced by skim::NueSkimmer::CopyMichelSlice(), and skim::NueSkimmer::CopyMichelTrack().

124  {
125  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
126  instanceName);
127  }
ModuleDescription moduleDescription_
Definition: EDProducer.h:115
template<typename PROD , BranchType B>
ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription,
std::string const &  instanceName 
) const
inherited

Definition at line 56 of file ProducerBase.h.

References art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
RVPInfo jmshower::RecoJMShower::GetRVPStats ( art::Handle< std::vector< rb::Cluster > >  slices,
unsigned int  ic 
)

Definition at line 5584 of file RecoJMShower_module.cc.

References std::abs(), rb::Cluster::Cell(), jmshower::RVPInfo::efrac_20slides, jmshower::RVPInfo::efrac_2sigmaHOR, jmshower::RVPInfo::efrac_2slides, jmshower::RVPInfo::efrac_4sigmaHOR, jmshower::RVPInfo::efrac_4slides, jmshower::RVPInfo::efrac_6slides, jmshower::RVPInfo::eiso_2sigmaHOR, jmshower::RVPInfo::eiso_4sigmaHOR, jmshower::RVPInfo::eiso_6sigmaHOR, rb::Cluster::ExtentPlane(), art::Ptr< T >::get(), rb::RecoHit::GeV(), MECModelEnuComparisons::i, rb::RecoHit::IsCalibrated(), geo::kX, calib::Calibrator::MakeRecoHit(), extractScale::mean, rb::Cluster::MeanXYZ(), rb::Cluster::MinPlane(), rb::Cluster::NCell(), nplanes, rb::Cluster::NXCell(), rb::Cluster::NYCell(), rb::CellHit::Plane(), rb::CellHit::View(), rb::Cluster::W(), rb::RecoHit::X(), rb::Cluster::XCell(), rb::RecoHit::Y(), and rb::Cluster::YCell().

Referenced by RecoShowers().

5585  {
5587 
5588  const rb::Cluster& slice = (*slices)[ic];
5589 
5590  RVPInfo rvp;
5591  rvp.efrac_2slides=1.;
5592  rvp.efrac_4slides=1.;
5593  rvp.efrac_6slides=1.;
5594  rvp.efrac_20slides=1.;
5595  rvp.efrac_2sigmaHOR=1.;
5596  rvp.efrac_4sigmaHOR=1.;
5597  rvp.eiso_2sigmaHOR=0.;
5598  rvp.eiso_4sigmaHOR=0.;
5599  rvp.eiso_6sigmaHOR=0.;
5600 
5601  const int lengthPlanes = slice.ExtentPlane();
5602  double recoEnergy=0.;
5603  for( unsigned int icell=0; icell<slice.NCell(); ++ icell){
5604  const art::Ptr<rb::CellHit>& chit = slice.Cell(icell);
5605  const TVector3 mean = slice.MeanXYZ();
5606  const double wx = (slice.NXCell() > 0) ? slice.W(slice.XCell(0).get()) : 0;
5607  const double wy = (slice.NYCell() > 0) ? slice.W(slice.YCell(0).get()) : 0;
5608 
5609  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, chit->View() == geo::kX ? wx : wy));
5610 
5611  if( !rhit.IsCalibrated() ) continue;
5612  recoEnergy += rhit.GeV();
5613  }//loop all cells
5614 
5615  //calculate per plane cell energy
5616  std::vector<double> PlaneSumEnergy;
5617  PlaneSumEnergy.clear();
5618  double totalE_road20=0.;
5619  double totalE_road40=0.;
5620  double totalE_road60=0.;
5621 
5622  for( int iplane=0; iplane<lengthPlanes; ++iplane ){
5623  double cellEnergy=0.;//sum of cells energy in each plane
5624  double cellX=0.;//energetic cell X in each plane
5625  double cellY=0.;//energetic cell Y in each plane
5626  double Ecell=0.;
5627  double EXcell=0.;
5628  double EYcell=0.;
5629 
5630  std::vector<double> Xxview;
5631  std::vector<double> Exview;
5632  std::vector<double> Yyview;
5633  std::vector<double> Eyview;
5634  Xxview.clear();
5635  Exview.clear();
5636  Yyview.clear();
5637  Eyview.clear();
5638 
5639  if( slice.NCell()>0 ){
5640  for( unsigned int icell=0; icell<slice.NCell(); ++ icell){
5641  const art::Ptr<rb::CellHit>& chit = slice.Cell(icell);
5642  const TVector3 mean = slice.MeanXYZ();
5643  const double wx = (slice.NXCell() > 0) ? slice.W(slice.XCell(0).get()) : 0;
5644  const double wy = (slice.NYCell() > 0) ? slice.W(slice.YCell(0).get()) : 0;
5645  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, chit->View() == geo::kX ? wx : wy));
5646 
5647  if( !rhit.IsCalibrated() ) continue;
5648 
5649  int Dplane=chit->Plane()-slice.MinPlane();
5650 
5651  if( Dplane==iplane ){
5652  cellEnergy += rhit.GeV();
5653  Ecell += rhit.GeV();
5654  EXcell += rhit.GeV()*rhit.X();
5655  EYcell += rhit.GeV()*rhit.Y();
5656 
5657  if( chit->View() == geo::kX ){
5658  Xxview.push_back(rhit.X());
5659  Exview.push_back(rhit.GeV());
5660  }//in X view
5661  else {
5662  Yyview.push_back(rhit.Y());
5663  Eyview.push_back(rhit.GeV());
5664  }//in Y view
5665  }
5666  }//loop cells
5667  }//has cells
5668  if( Ecell>0. ){
5669  cellX = EXcell/Ecell;
5670  cellY = EYcell/Ecell;
5671  }
5672 
5673  PlaneSumEnergy.push_back(cellEnergy);
5674 
5675  int Nxviews=Xxview.size();
5676  int Nyviews=Yyview.size();
5677  if( Nxviews>0 ){
5678  for( int ixv=0; ixv<Nxviews; ++ixv ){
5679  if( std::abs(Xxview[ixv]-cellX) < 2.*2. ) totalE_road20+=Exview[ixv];
5680  if( std::abs(Xxview[ixv]-cellX) < 2.*4. ) totalE_road40+=Exview[ixv];
5681  if( std::abs(Xxview[ixv]-cellX) < 2.*6. ) totalE_road60+=Exview[ixv];
5682  }
5683  }//in X view
5684  if( Nyviews>0 ){
5685  for( int iyv=0; iyv<Nyviews; ++iyv ){
5686  if( std::abs(Yyview[iyv]-cellY) < 2.*2. ) totalE_road20+=Eyview[iyv];
5687  if( std::abs(Yyview[iyv]-cellY) < 2.*4. ) totalE_road40+=Eyview[iyv];
5688  if( std::abs(Yyview[iyv]-cellY) < 2.*6. ) totalE_road60+=Eyview[iyv];
5689  }
5690  }//in Y view
5691 
5692 /*
5693  if( slice.NCell()>0 ){
5694  for( unsigned int icell=0; icell<slice.NCell(); ++ icell){
5695  const art::Ptr<rb::CellHit>& chit = slice.Cell(icell);
5696  const TVector3 mean = slice.MeanXYZ();
5697  const double wx = (slice.NXCell() > 0) ? slice.W(slice.XCell(0).get()) : 0;
5698  const double wy = (slice.NYCell() > 0) ? slice.W(slice.YCell(0).get()) : 0;
5699  const rb::RecoHit rhit = cal->MakeRecoHit(*chit,
5700  chit->View() == geo::kX ? wx : wy);
5701  if( !rhit.IsCalibrated() ) continue;
5702  int Dplane=chit->Plane()-slice.MinPlane();
5703  if( Dplane==iplane ){
5704  double dx=std::abs(rhit.X()-cellX);
5705  double dy=std::abs(rhit.Y()-cellY);
5706  if( chit->View()==geo::kX ){
5707  if( dx<2.*2. ) totalE_road20+=rhit.GeV();
5708  if( dx<2.*4. ) totalE_road40+=rhit.GeV();
5709  }//in X view
5710  else{
5711  if( dy<2.*2. ) totalE_road20+=rhit.GeV();
5712  if( dy<2.*4. ) totalE_road40+=rhit.GeV();
5713  }//in Y view
5714  }
5715  }//loop all cells
5716  }//has cells
5717  */
5718 
5719  }//loop planes
5720 
5721  double efrac_plane2=-1.;
5722  double efrac_plane4=-1.;
5723  double efrac_plane6=-1.;
5724  double Eplane20=0.;
5725  int nplanes=PlaneSumEnergy.size();
5726  if (recoEnergy > 0.0){
5727  for( int i=0; i<nplanes; ++i ){
5728  if( i<20 ) Eplane20 += PlaneSumEnergy[i];
5729  //maximal energy fraction in 2 sliding window
5730  if( i<nplanes-2 ){
5731  double e_frac2=(PlaneSumEnergy[i]+PlaneSumEnergy[i+1])/recoEnergy;
5732  if( efrac_plane2<e_frac2 ) efrac_plane2=e_frac2;
5733  }
5734  else if( nplanes<=2 ) efrac_plane2=1.;
5735  //maximal energy fraction in 4 sliding window
5736  if( i<nplanes-4 ){
5737  double e_frac4=(PlaneSumEnergy[i]+PlaneSumEnergy[i+1]+PlaneSumEnergy[i+2]+PlaneSumEnergy[i+3])/recoEnergy;
5738  if( efrac_plane4<e_frac4 ) efrac_plane4=e_frac4;
5739  }
5740  else if( nplanes<=4 ) efrac_plane4=1.;
5741  //maximal energy fraction in 6 sliding window
5742  if( i<nplanes-6 ){
5743  double e_frac6=(PlaneSumEnergy[i]+PlaneSumEnergy[i+1]+PlaneSumEnergy[i+2]+PlaneSumEnergy[i+3]+PlaneSumEnergy[i+4]+PlaneSumEnergy[i+5])/recoEnergy;
5744  if( efrac_plane6<e_frac6 ) efrac_plane6=e_frac6;
5745  }
5746  else if( nplanes<=6 ) efrac_plane6=1.;
5747  }
5748  }
5749  if( recoEnergy>0.0 ){
5750  rvp.efrac_2sigmaHOR=totalE_road20/recoEnergy;
5751  rvp.efrac_4sigmaHOR=totalE_road40/recoEnergy;
5752  rvp.efrac_20slides=Eplane20/recoEnergy;
5753  rvp.efrac_2slides=efrac_plane2;
5754  rvp.efrac_4slides=efrac_plane4;
5755  rvp.efrac_6slides=efrac_plane6;
5756  }
5757  else{
5758  rvp.efrac_2sigmaHOR=-1.0;
5759  rvp.efrac_4sigmaHOR=-1.0;
5760  rvp.efrac_20slides=-1.0;
5761  rvp.efrac_2slides=efrac_plane2;
5762  rvp.efrac_4slides=efrac_plane4;
5763  rvp.efrac_6slides=efrac_plane6;
5764  }
5765 
5766  if( totalE_road20>0. ) rvp.eiso_2sigmaHOR=(recoEnergy-totalE_road20)/totalE_road20;
5767  if( totalE_road40>0. ) rvp.eiso_4sigmaHOR=(recoEnergy-totalE_road40)/totalE_road40;
5768  if( totalE_road60>0. ) rvp.eiso_6sigmaHOR=(recoEnergy-totalE_road60)/totalE_road60;
5769  return rvp;
5770  }
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Cluster.cxx:121
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
Vertical planes which measure X.
Definition: PlaneGeo.h:28
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
A collection of associated CellHits.
Definition: Cluster.h:47
TVector3 MeanXYZ(rb::AveragingScheme=kDefaultScheme) const
Definition: Cluster.cxx:538
float abs(float number)
Definition: d0nt_math.hpp:39
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
int nplanes
Definition: geom.C:145
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
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
T const * get() const
Definition: Ptr.h:321
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
Definition: FillPIDs.h:15
TVector3 jmshower::RecoJMShower::GetShwStop ( rb::Prong  shw)

Definition at line 3006 of file RecoJMShower_module.cc.

References rb::Prong::Dir(), rb::Prong::Start(), and rb::Prong::TotalLength().

Referenced by Energy(), GetDedxInvLongLL(), GetDedxLongLL(), GetDedxTransLL(), GetTrkHitPath(), GetTrkHitPos(), GetTrkPlaneCell(), GetTrkPlaneDistToEdge(), GetTrkPlanePos(), Radius(), and RecoShowers().

3006  {
3007  TVector3 shwdir=shw.Dir();
3008  TVector3 shwstart=shw.Start();
3009  TVector3 shwend = shwstart + shw.TotalLength()*shwdir;
3010  return shwend;
3011  }
virtual TVector3 Start() const
Definition: Prong.h:73
virtual double TotalLength() const
Distance along prong to reach last cell hit.
Definition: Prong.cxx:186
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
TVector3 jmshower::RecoJMShower::GetShwStop ( art::Ptr< rb::Prong shw)

Definition at line 3013 of file RecoJMShower_module.cc.

References rb::Prong::Dir(), rb::Prong::Start(), and rb::Prong::TotalLength().

3013  {
3014  TVector3 shwdir=shw->Dir();
3015  TVector3 shwstart=shw->Start();
3016  TVector3 shwend = shwstart + shw->TotalLength()*shwdir;
3017  return shwend;
3018  }
virtual TVector3 Start() const
Definition: Prong.h:73
virtual double TotalLength() const
Distance along prong to reach last cell hit.
Definition: Prong.cxx:186
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
unsigned int jmshower::RecoJMShower::GetTrkCPlaneCell ( unsigned int  plane,
unsigned int  i 
)

Definition at line 3286 of file RecoJMShower_module.cc.

References rb::CellHit::Cell(), cellPos(), PandAna.Demos.tute_pid_validation::cid, fHitEWeightMap, fShowerPClusterCol, geom(), MECModelEnuComparisons::i, makeTrainCVSamples::int, findDuplicateFiles::key, geo::PlaneGeo::Ncells(), rb::CellHit::Plane(), geo::GeometryBase::Plane(), and NDAPDHVSetting::plane.

Referenced by RecoShowers().

3286  {//calculate centroid of a plane
3288  double sumid = 0;
3289  double gev = 0;
3290  for (unsigned int nh=0;nh<fShowerPClusterCol[i][plane].NCell();nh++) {
3291  art::Ptr<rb::CellHit> cellhit = fShowerPClusterCol[i][plane].Cell(nh);
3292  geo::OfflineChan key(cellhit->Plane(),cellhit->Cell());
3293  double cellEnergy = fHitEWeightMap[key][i];
3294  double cellPos = (double)cellhit->Cell() + 0.5;
3295  sumid += cellEnergy*cellPos;
3296  gev += cellEnergy;
3297  }
3298  if(gev<0.00001)return 9999;
3299  double avgid = sumid/gev;
3300  int cid = (int)avgid;
3301  if(cid<0)cid=0;
3302  if(cid>(int)geom->Plane(plane)->Ncells()-1)cid=(int)geom->Plane(plane)->Ncells()-1;
3303  return (unsigned int)cid;
3304  }
std::map< geo::OfflineChan, std::vector< double > > fHitEWeightMap
unsigned short Plane() const
Definition: CellHit.h:39
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
const PlaneGeo * Plane(unsigned int i) const
std::vector< std::vector< jmshower::PCluster > > fShowerPClusterCol
unsigned short Cell() const
Definition: CellHit.h:40
void cellPos()
Definition: cellShifts.C:26
A (plane, cell) pair.
Definition: OfflineChan.h:17
void geom(int which=0)
Definition: geom.C:163
double jmshower::RecoJMShower::GetTrkDistToEdge ( rb::Prong  shower)

Definition at line 3451 of file RecoJMShower_module.cc.

References std::abs(), fShowerPClusterCol, GetTrkPlaneDistToEdge(), MECModelEnuComparisons::i, rb::Cluster::ID(), and NDAPDHVSetting::plane.

Referenced by RecoShowers().

3451  {
3452  int i = shower.ID();
3453  double mindist = 999999.;
3454  for(unsigned int plane=0;plane<fShowerPClusterCol[i].size();plane++){
3455  TVector3 disttoedge = RecoJMShower::GetTrkPlaneDistToEdge(fShowerPClusterCol[i][plane].Plane(), shower);
3456  double minxyz=999999.;
3457  for(int ii=0;ii<3;ii++){
3458  if(std::abs(disttoedge[ii])<minxyz)minxyz=std::abs(disttoedge[ii]);
3459  }
3460  if(minxyz<mindist)mindist=minxyz;
3461  }
3462  return mindist;
3463  }
std::vector< std::vector< jmshower::PCluster > > fShowerPClusterCol
float abs(float number)
Definition: d0nt_math.hpp:39
const int ID() const
Definition: Cluster.h:75
TVector3 GetTrkPlaneDistToEdge(unsigned int plane, art::Ptr< rb::Prong > trk)
double jmshower::RecoJMShower::GetTrkHitPath ( unsigned int  plane,
unsigned int  cell,
rb::Prong  shower 
)

Definition at line 3466 of file RecoJMShower_module.cc.

References std::abs(), geo::PlaneGeo::Cell(), om::cout, rb::Prong::Dir(), e, allTimeWatchdog::endl, fShowerPClusterCol, geom(), geo::CellGeo::GetCenter(), GetShwStop(), geo::CellGeo::HalfD(), geo::CellGeo::HalfL(), MECModelEnuComparisons::i, rb::Cluster::ID(), Mag(), geo::GeometryBase::Plane(), gen_flatrecord::size, std::sqrt(), and rb::Prong::Start().

Referenced by GetCellTransDedx(), and GetPlaneDedx().

3466  {
3467  int i = shower.ID();
3468  Double_t cellXYZ[3]={0,0,0}; // xyz[0] = x, xyz[1] = y, xyz[2] = z
3470  geom->Plane(plane)->Cell(cell)->GetCenter(cellXYZ);
3471  TVector3 trkdir=shower.Dir();
3472  TVector3 trkstart=shower.Start();
3473  TVector3 trkstop=RecoJMShower::GetShwStop(shower);
3474  double trklength = std::abs((trkstart - trkstop).Mag());
3475  double cellD=2*geom->Plane(plane)->Cell(cell)->HalfD();
3476  double cellL=2*geom->Plane(plane)->Cell(cell)->HalfL();
3477  if(cellL<cellD||cellL<0.0001){
3478  std::cout<<"Geometry error! "<<std::endl;
3479  return 9999;
3480  }
3481  if(std::abs(trkdir.z())< 1e-6||fShowerPClusterCol[i].size()==1){
3482 // std::cout<<"Small Z Dir! "<<std::endl;
3483 // std::cout<<"N Plane = "<<shower.ExtentPlane()<<std::endl;
3484 // std::cout<<"Path = "<<trklength<<std::endl;
3485  return trklength;
3486  }
3487  double trkhitpath=cellD*(sqrt(trkdir.x()*trkdir.x()+trkdir.y()*trkdir.y()+trkdir.z()*trkdir.z())/std::abs(trkdir.z()));
3488  return trkhitpath;
3489  }
double HalfL() const
Definition: CellGeo.cxx:198
double HalfD() const
Definition: CellGeo.cxx:205
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 GetShwStop(rb::Prong shw)
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
T sqrt(T number)
Definition: d0nt_math.hpp:156
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
std::vector< std::vector< jmshower::PCluster > > fShowerPClusterCol
float abs(float number)
Definition: d0nt_math.hpp:39
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
OStream cout
Definition: OStream.cxx:6
void geom(int which=0)
Definition: geom.C:163
float Mag() const
const int ID() const
Definition: Cluster.h:75
Float_t e
Definition: plot.C:35
TVector3 jmshower::RecoJMShower::GetTrkHitPos ( unsigned int  plane,
unsigned int  cell,
art::Ptr< rb::Prong trk 
)

Definition at line 3022 of file RecoJMShower_module.cc.

References std::abs(), geo::PlaneGeo::Cell(), rb::Prong::Dir(), geo::GeometryBase::DistToElectronics(), geom(), geo::CellGeo::GetCenter(), GetShwStop(), geo::kX, geo::kY, geo::GeometryBase::Plane(), rb::Prong::Start(), POTSpillRate::view, and geo::PlaneGeo::View().

Referenced by RecoShowers().

3022  {
3023  Double_t cellXYZ[3]={0,0,0}; //declare container to hold cell coordinate
3025  geom->Plane(plane)->Cell(cell)->GetCenter(cellXYZ); //get cell center
3026 
3027  //get track direction, start, end
3028  TVector3 trkdir=trk->Dir();
3029  TVector3 trkstart=trk->Start();
3030  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3031 // double cellL=2*geom->Plane(plane)->Cell(cell)->HalfL();
3032  geo::View_t view = geom->Plane(plane)->View();
3033  double celldist = 0;
3034  double trkx = 0;
3035  double trky = 0;
3036  //calculate the tracks x position at the cell's z coordinate
3037  trkx = trkstart.x()+(trkdir.x()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3038  //calculate the tracks y position at the cell's z coordinate
3039  trky = trkstart.y()+(trkdir.y()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3040 // double trkdtr0 = 0;
3041  double trkdtr = 0;
3042  if(view == geo::kX){
3043  //calculate distance from cell x coordinate to track x coordninate at same z position
3044  //NOTE: Perpendicular distance to track may be a better metric as this will be incorrect for steep tracks
3045  celldist=std::abs((cellXYZ[0]-trkstart.x())-((cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z())));
3046 // trkdtr0 = geom->Plane(plane)->Cell(cell)->DistToReadOut(trky); //distance to readout for cell
3047  trkdtr = geom->DistToElectronics(trky, *geom->Plane(plane)->Cell(cell));
3048 
3049  }else if(view == geo::kY){
3050  //calculate distance from cell x coordinate to track x coordninate at same z position
3051  celldist=std::abs((cellXYZ[1]-trkstart.y())-((cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z())));
3052 // trkdtr0 = geom->Plane(plane)->Cell(cell)->DistToReadOut(trkx); //distance to readout for cell
3053  trkdtr = geom->DistToElectronics(trkx, *geom->Plane(plane)->Cell(cell));
3054  }
3055  if(trkdtr<0)trkdtr=0;
3056  //return results
3057  TVector3 def(-1000,-1000,-1000);
3058  TVector3 trkhitposxcell(trkdtr,trky,celldist);
3059  TVector3 trkhitposycell(trkdtr,trkx,celldist);
3060  if(view == geo::kX){return trkhitposxcell;}
3061  if(view == geo::kY){return trkhitposycell;}
3062  return def;
3063  }
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 GetShwStop(rb::Prong shw)
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
virtual TVector3 Start() const
Definition: Prong.h:73
float abs(float number)
Definition: d0nt_math.hpp:39
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
double DistToElectronics(double localz, const CellGeo &cell) const
get distance from local z position in cell to apd in cm, including pigtail
void geom(int which=0)
Definition: geom.C:163
TVector3 jmshower::RecoJMShower::GetTrkHitPos ( unsigned int  plane,
unsigned int  cell,
rb::Prong  trk 
)

Definition at line 3066 of file RecoJMShower_module.cc.

References std::abs(), geo::PlaneGeo::Cell(), rb::Prong::Dir(), geo::GeometryBase::DistToElectronics(), geom(), geo::CellGeo::GetCenter(), GetShwStop(), geo::kX, geo::kY, geo::GeometryBase::Plane(), rb::Prong::Start(), POTSpillRate::view, and geo::PlaneGeo::View().

3066  {
3067  Double_t cellXYZ[3]={0,0,0};
3069  geom->Plane(plane)->Cell(cell)->GetCenter(cellXYZ);
3070  TVector3 trkdir=trk.Dir();
3071  TVector3 trkstart=trk.Start();
3072  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3073 // double cellL=2*geom->Plane(plane)->Cell(cell)->HalfL();
3074  geo::View_t view = geom->Plane(plane)->View();
3075  double celldist = 0;
3076  double trkx = 0;
3077  double trky = 0;
3078  trkx = trkstart.x()+(trkdir.x()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3079  trky = trkstart.y()+(trkdir.y()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3080 
3081 // double trkdtr0 = 0;
3082  double trkdtr = 0;
3083  if(view == geo::kX){
3084  celldist=std::abs((cellXYZ[0]-trkstart.x())-((cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z())));
3085 // trkdtr = geom->Plane(plane)->Cell(cell)->DistToReadOut(trky);
3086  trkdtr = geom->DistToElectronics(trky, *geom->Plane(plane)->Cell(cell));
3087  }else if(view == geo::kY){
3088  celldist=std::abs((cellXYZ[1]-trkstart.y())-((cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z())));
3089 // trkdtr = geom->Plane(plane)->Cell(cell)->DistToReadOut(trkx);
3090  trkdtr = geom->DistToElectronics(trkx, *geom->Plane(plane)->Cell(cell));
3091  }
3092  if(trkdtr<0)trkdtr=0;
3093  TVector3 def(-1000,-1000,-1000);
3094  TVector3 trkhitposxcell(trkdtr,trky,celldist);
3095  TVector3 trkhitposycell(trkdtr,trkx,celldist);
3096  if(view == geo::kX){return trkhitposxcell;}
3097  if(view == geo::kY){return trkhitposycell;}
3098  return def;
3099  }
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 GetShwStop(rb::Prong shw)
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
virtual TVector3 Start() const
Definition: Prong.h:73
float abs(float number)
Definition: d0nt_math.hpp:39
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
double DistToElectronics(double localz, const CellGeo &cell) const
get distance from local z position in cell to apd in cm, including pigtail
void geom(int which=0)
Definition: geom.C:163
TVector3 jmshower::RecoJMShower::GetTrkHitPos ( unsigned int  plane,
unsigned int  cell,
TVector3  trkdir,
TVector3  trkstart,
TVector3  trkstop 
)

Definition at line 3101 of file RecoJMShower_module.cc.

References std::abs(), geo::PlaneGeo::Cell(), geo::GeometryBase::DistToElectronics(), geom(), geo::CellGeo::GetCenter(), geo::kX, geo::kY, geo::GeometryBase::Plane(), POTSpillRate::view, and geo::PlaneGeo::View().

3101  {
3102  Double_t cellXYZ[3]={0,0,0}; // XYZ[0] = x, XYZ[1] = y, XYZ[2] = z
3104  geom->Plane(plane)->Cell(cell)->GetCenter(cellXYZ);
3105 // double cellL=2*geom->Plane(plane)->Cell(cell)->HalfL();
3106  geo::View_t view = geom->Plane(plane)->View();
3107  double celldist = 0;
3108  double trkx = 0;
3109  double trky = 0;
3110  trkx = trkstart.x()+(trkdir.x()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3111  trky = trkstart.y()+(trkdir.y()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3112 // double trkdtr0 = 0;
3113  double trkdtr = 0;
3114  if(view == geo::kX){
3115  celldist=std::abs((cellXYZ[0]-trkstart.x())-((cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z())));
3116 // trkdtr0 = geom->Plane(plane)->Cell(cell)->DistToReadOut(trky);
3117  trkdtr = geom->DistToElectronics(trky, *geom->Plane(plane)->Cell(cell));
3118  }else if(view == geo::kY){
3119  celldist=std::abs((cellXYZ[1]-trkstart.y())-((cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z())));
3120 // trkdtr0 = geom->Plane(plane)->Cell(cell)->DistToReadOut(trkx);
3121  trkdtr = geom->DistToElectronics(trkx, *geom->Plane(plane)->Cell(cell));
3122  }
3123  if(trkdtr<0)trkdtr=0;
3124  TVector3 def(-1000,-1000,-1000);
3125  TVector3 trkhitposxcell(trkdtr,trky,celldist);
3126  TVector3 trkhitposycell(trkdtr,trkx,celldist);
3127  if(view == geo::kX){return trkhitposxcell;}
3128  if(view == geo::kY){return trkhitposycell;}
3129  return def;
3130  }
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
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
double DistToElectronics(double localz, const CellGeo &cell) const
get distance from local z position in cell to apd in cm, including pigtail
void geom(int which=0)
Definition: geom.C:163
unsigned int jmshower::RecoJMShower::GetTrkPlaneCell ( unsigned int  plane,
art::Ptr< rb::Prong trk 
)

Definition at line 3177 of file RecoJMShower_module.cc.

References std::abs(), geo::PlaneGeo::Cell(), rb::Prong::Dir(), geom(), geo::CellGeo::GetCenter(), GetShwStop(), geo::CellGeo::HalfW(), makeTrainCVSamples::int, geo::kX, geo::kY, geo::PlaneGeo::Ncells(), geo::GeometryBase::Plane(), rb::Prong::Start(), POTSpillRate::view, and geo::PlaneGeo::View().

3177  {
3178 
3179  Double_t cellXYZ[3]={0,0,0};
3180  Double_t cellXYZ1[3];
3181  Double_t cellXYZmax[3];
3183  geom->Plane(plane)->Cell(0)->GetCenter(cellXYZ);
3184  geom->Plane(plane)->Cell(geom->Plane(plane)->Ncells()-1)->GetCenter(cellXYZ1);
3185  geom->Plane(plane)->Cell(geom->Plane(plane)->Ncells()-1)->GetCenter(cellXYZmax);
3186  double cellW = 2*geom->Plane(plane)->Cell(0)->HalfW();
3187  double cellwx = (cellXYZ1[0] - cellXYZ[0] + cellW)/geom->Plane(plane)->Ncells();
3188  double cellwy = (cellXYZ1[1] - cellXYZ[1] + cellW)/geom->Plane(plane)->Ncells();
3189  TVector3 trkdir=trk->Dir();
3190  TVector3 trkstart=trk->Start();
3191  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3192  double trkplanex=99999,trkplaney=99999;//, trkplanez=0;
3193 // double trkplanex0=99999,trkplaney0=99999;
3194 // double trkplanex1=99999,trkplaney1=99999;
3195  if(std::abs(trkdir.z())>0.0001){
3196  trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3197  trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3198 // trkplanex1 = cellXYZmax[0];
3199 // trkplaney1 = cellXYZmax[1];
3200  }else{
3201  trkplanex = (trkstart.x()+ trkstop.x())/2.0;
3202  trkplaney = (trkstart.y()+ trkstop.y())/2.0;
3203  }
3204 // trkplanez = cellXYZ[2];
3205  geo::View_t view = geom->Plane(plane)->View();
3206 /*
3207  if(view == geo::kX){
3208  std::cout<<"x view plane"<<std::endl;
3209  std::cout<<"x1-x0,cellW "<<" "<<cellXYZ1[0] - cellXYZ[0]<<" "<<cellW<<std::endl;
3210  }else if(view == geo::kY){
3211  std::cout<<"y view plane"<<std::endl;
3212  std::cout<<"y1-y0,cellW "<<" "<<cellXYZ1[1] - cellXYZ[1]<<" "<<cellW<<std::endl;
3213  }
3214 */
3215  int trkcellId = -9999;
3216 // int trkcellId0 = -9999;
3217 // int trkcellId1 = -9999;
3218 // int trkcellIdmax = geom->Plane(plane)->Ncells()-1;
3219  if(view == geo::kX)trkcellId = int((trkplanex+cellwx/2-cellXYZ[0])/cellwx);
3220  if(view == geo::kY)trkcellId = int((trkplaney+cellwy/2-cellXYZ[1])/cellwy);
3221 // if(view == geo::kX)trkcellId0 = int((trkplanex0+cellwx/2-cellXYZ[0])/cellwx);
3222 // if(view == geo::kY)trkcellId0 = int((trkplaney0+cellwy/2-cellXYZ[1])/cellwy);
3223 // if(view == geo::kX)trkcellId1 = int((trkplanex1+cellwx/2-cellXYZ[0])/cellwx);
3224 // if(view == geo::kY)trkcellId1 = int((trkplaney1+cellwy/2-cellXYZ[1])/cellwy);
3225  return unsigned(trkcellId);
3226  }
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 GetShwStop(rb::Prong shw)
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
double HalfW() const
Definition: CellGeo.cxx:191
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
float abs(float number)
Definition: d0nt_math.hpp:39
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
void geom(int which=0)
Definition: geom.C:163
unsigned int jmshower::RecoJMShower::GetTrkPlaneCell ( unsigned int  plane,
rb::Prong  trk 
)

Definition at line 3228 of file RecoJMShower_module.cc.

References std::abs(), geo::PlaneGeo::Cell(), rb::Prong::Dir(), geom(), geo::CellGeo::GetCenter(), GetShwStop(), geo::CellGeo::HalfW(), makeTrainCVSamples::int, geo::kX, geo::kY, geo::PlaneGeo::Ncells(), geo::GeometryBase::Plane(), rb::Prong::Start(), POTSpillRate::view, and geo::PlaneGeo::View().

3228  {
3229 
3230  Double_t cellXYZ[3]={0,0,0};
3231  Double_t cellXYZ1[3];
3232  Double_t cellXYZmax[3];
3234  geom->Plane(plane)->Cell(0)->GetCenter(cellXYZ);
3235  geom->Plane(plane)->Cell(geom->Plane(plane)->Ncells()-1)->GetCenter(cellXYZ1);
3236  geom->Plane(plane)->Cell(geom->Plane(plane)->Ncells()-1)->GetCenter(cellXYZmax);
3237  double cellW = 2*geom->Plane(plane)->Cell(0)->HalfW();
3238  double cellwx = (cellXYZ1[0] - cellXYZ[0] + cellW)/geom->Plane(plane)->Ncells();
3239  double cellwy = (cellXYZ1[1] - cellXYZ[1] + cellW)/geom->Plane(plane)->Ncells();
3240  TVector3 trkdir=trk.Dir();
3241  TVector3 trkstart=trk.Start();
3242  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3243  double trkplanex=99999,trkplaney=99999;//, trkplanez=0;
3244 // double trkplanex0=99999,trkplaney0=99999;
3245 // double trkplanex1=99999,trkplaney1=99999;
3246  if(std::abs(trkdir.z())>0.0001){
3247  trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3248  trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3249 // trkplanex0 = cellXYZ[0];
3250 // trkplaney0 = cellXYZ[1];
3251 // trkplanex1 = cellXYZmax[0];
3252 // trkplaney1 = cellXYZmax[1];
3253  }else{
3254  trkplanex = (trkstart.x()+trkstop.x())/2.0;
3255  trkplaney = (trkstart.y()+trkstop.y())/2.0;
3256  }
3257 // trkplanez = cellXYZ[2];
3258  geo::View_t view = geom->Plane(plane)->View();
3259 
3260 /*
3261  if(view == geo::kX){
3262  std::cout<<"x view plane"<<std::endl;
3263  std::cout<<"x1-x0,cellW "<<" "<<cellXYZ1[0] - cellXYZ[0]<<" "<<cellW<<std::endl;
3264  }else if(view == geo::kY){
3265  std::cout<<"y view plane"<<std::endl;
3266  std::cout<<"y1-y0,cellW "<<" "<<cellXYZ1[1] - cellXYZ[1]<<" "<<cellW<<std::endl;
3267  }
3268 */
3269 
3270  int trkcellId = -9999;
3271 // int trkcellId0 = -9999;
3272 // int trkcellId1 = -9999;
3273 // int trkcellIdmax = geom->Plane(plane)->Ncells()-1;
3274 
3275  if(view == geo::kX)trkcellId = int((trkplanex+cellwx/2-cellXYZ[0])/cellwx);
3276  if(view == geo::kY)trkcellId = int((trkplaney+cellwy/2-cellXYZ[1])/cellwy);
3277 
3278 // if(view == geo::kX)trkcellId0 = int((trkplanex0+cellwx/2-cellXYZ[0])/cellwx);
3279 // if(view == geo::kY)trkcellId0 = int((trkplaney0+cellwy/2-cellXYZ[1])/cellwy);
3280 //
3281 // if(view == geo::kX)trkcellId1 = int((trkplanex1+cellwx/2-cellXYZ[0])/cellwx);
3282 // if(view == geo::kY)trkcellId1 = int((trkplaney1+cellwy/2-cellXYZ[1])/cellwy);
3283  return unsigned(trkcellId);
3284  }
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 GetShwStop(rb::Prong shw)
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
double HalfW() const
Definition: CellGeo.cxx:191
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
float abs(float number)
Definition: d0nt_math.hpp:39
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
void geom(int which=0)
Definition: geom.C:163
TVector3 jmshower::RecoJMShower::GetTrkPlaneDistToEdge ( unsigned int  plane,
art::Ptr< rb::Prong trk 
)

Definition at line 3306 of file RecoJMShower_module.cc.

References std::abs(), geo::PlaneGeo::Cell(), rb::Prong::Dir(), geom(), geo::CellGeo::GetCenter(), GetShwStop(), geo::GeometryBase::Plane(), pXEdge1, pXEdge2, pXEdge3, pXEdge4, pYEdge1, pYEdge2, pYEdge3, pYEdge4, pZEdge1, pZEdge2, pZEdge3, pZEdge4, and rb::Prong::Start().

Referenced by GetTrkDistToEdge(), and IsFiducial().

3306  {
3307  Double_t cellXYZ[3]={0,0,0};
3309  geom->Plane(plane)->Cell(0)->GetCenter(cellXYZ);
3310  TVector3 trkdir=trk->Dir();
3311  TVector3 trkstart=trk->Start();
3312  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3313  double trkplanex=99999,trkplaney=99999, trkplanez=0;
3314  if(std::abs(trkdir.z())>0.0001){
3315  trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3316  trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3317  }
3318  trkplanez = cellXYZ[2];
3319 
3320  double xe[4]={9999,9999,9999,9999};
3321  double ye[4]={9999,9999,9999,9999};
3322  double ze[4]={9999,9999,9999,9999};
3323 
3324  xe[0] = std::abs(pXEdge1-trkplanex);
3325  xe[1] = std::abs(pXEdge2-trkplanex);
3326  xe[2] = std::abs(pXEdge3-trkplanex);
3327  xe[3] = std::abs(pXEdge4-trkplanex);
3328 
3329  ye[0] = std::abs(pYEdge1-trkplaney);
3330  ye[1] = std::abs(pYEdge2-trkplaney);
3331  ye[2] = std::abs(pYEdge3-trkplaney);
3332  ye[3] = std::abs(pYEdge4-trkplaney);
3333 
3334  ze[0] = std::abs(pZEdge1-trkplanez);
3335  ze[1] = std::abs(pZEdge2-trkplanez);
3336  ze[2] = std::abs(pZEdge3-trkplanez);
3337  ze[3] = std::abs(pZEdge4-trkplanez);
3338 
3339  double xtoedge = 9999.;
3340  double ytoedge = 9999.;
3341  double ztoedge = 9999.;
3342 
3343  for(int ii=0;ii<4;ii++){
3344  if(xtoedge>xe[ii])xtoedge=xe[ii];
3345  if(ytoedge>ye[ii])ytoedge=ye[ii];
3346  if(ztoedge>ze[ii])ztoedge=ze[ii];
3347  }
3348 
3349 // double xtoedge0 = std::abs(geom->DetHalfWidth()-std::abs(trkplanex));
3350 // double ytoedge0 = std::abs(geom->DetHalfHeight()-std::abs(trkplaney));
3351 // double ztoedge0 = std::min(std::abs(geom->DetLength()-std::abs(trkplanez)),trkplanez);
3352 // std::cout<<"Det W/2,H/2,L "<<geom->DetHalfWidth()<<" "<<geom->DetHalfHeight()<<" "<<geom->DetLength()<<std::endl;
3353 // double cellL = 2*geom->Plane(plane)->Cell(0)->HalfL();
3354 // double cellW = 2*geom->Plane(plane)->Cell(0)->HalfW();
3355 // double cellD = 2*geom->Plane(plane)->Cell(0)->HalfD();
3356 // geo::View_t view = geom->Plane(plane)->View();
3357 /*
3358  if(view == geo::kX){
3359  std::cout<<"x view plane"<<std::endl;
3360  }else if(view == geo::kY){
3361  std::cout<<"y view plane"<<std::endl;
3362  }
3363  std::cout<<"cell L, cellW, cell D "<<cellL<<" "<<cellW<<" "<<cellD<<std::endl;
3364 */
3365 // std::cout<<"To Edge X,Y,Z "<<xtoedge<<" "<<ytoedge<<" "<<ztoedge<<std::endl;
3366  TVector3 toedge(xtoedge, ytoedge, ztoedge);
3367  return toedge;
3368  }
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 GetShwStop(rb::Prong shw)
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
float abs(float number)
Definition: d0nt_math.hpp:39
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
void geom(int which=0)
Definition: geom.C:163
TVector3 jmshower::RecoJMShower::GetTrkPlaneDistToEdge ( unsigned int  plane,
rb::Prong  trk 
)

Definition at line 3371 of file RecoJMShower_module.cc.

References std::abs(), geo::PlaneGeo::Cell(), rb::Prong::Dir(), geom(), geo::CellGeo::GetCenter(), GetShwStop(), geo::GeometryBase::Plane(), pXEdge1, pXEdge2, pXEdge3, pXEdge4, pYEdge1, pYEdge2, pYEdge3, pYEdge4, pZEdge1, pZEdge2, pZEdge3, pZEdge4, and rb::Prong::Start().

3371  {
3372  Double_t cellXYZ[3]={0,0,0};
3374  geom->Plane(plane)->Cell(0)->GetCenter(cellXYZ);
3375  TVector3 trkdir=trk.Dir();
3376  TVector3 trkstart=trk.Start();
3377  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3378  double trkplanex=99999,trkplaney=99999, trkplanez=0;
3379  if(std::abs(trkdir.z())>0.0001){
3380  trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3381  trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3382  }
3383  trkplanez = cellXYZ[2];
3384 
3385  double xe[4]={9999,9999,9999,9999};
3386  double ye[4]={9999,9999,9999,9999};
3387  double ze[4]={9999,9999,9999,9999};
3388 
3389  xe[0] = std::abs(pXEdge1-trkplanex);
3390  xe[1] = std::abs(pXEdge2-trkplanex);
3391  xe[2] = std::abs(pXEdge3-trkplanex);
3392  xe[3] = std::abs(pXEdge4-trkplanex);
3393 
3394  ye[0] = std::abs(pYEdge1-trkplaney);
3395  ye[1] = std::abs(pYEdge2-trkplaney);
3396  ye[2] = std::abs(pYEdge3-trkplaney);
3397  ye[3] = std::abs(pYEdge4-trkplaney);
3398 
3399  ze[0] = std::abs(pZEdge1-trkplanez);
3400  ze[1] = std::abs(pZEdge2-trkplanez);
3401  ze[2] = std::abs(pZEdge3-trkplanez);
3402  ze[3] = std::abs(pZEdge4-trkplanez);
3403 
3404  double xtoedge = 9999.;
3405  double ytoedge = 9999.;
3406  double ztoedge = 9999.;
3407 
3408  for(int ii=0;ii<4;ii++){
3409  if(xtoedge>xe[ii])xtoedge=xe[ii];
3410  if(ytoedge>ye[ii])ytoedge=ye[ii];
3411  if(ztoedge>ze[ii])ztoedge=ze[ii];
3412  }
3413 
3414 // double xtoedge0 = std::abs(geom->DetHalfWidth()-std::abs(trkplanex));
3415 // double ytoedge0 = std::abs(geom->DetHalfHeight()-std::abs(trkplaney));
3416 // double ztoedge0 = std::min(std::abs(geom->DetLength()-std::abs(trkplanez)),trkplanez);
3417 
3418 // double xtoedge = geom->DetHalfWidth()-std::abs(trkplanex);
3419 // double ytoedge = geom->DetHalfHeight()-std::abs(trkplaney);
3420 // double ztoedge = geom->DetLength()-std::abs(trkplanez);
3421 // std::cout<<"Det W/2,H/2,L "<<geom->DetHalfWidth()<<" "<<geom->DetHalfHeight()<<" "<<geom->DetLength()<<std::endl;
3422 // double cellL = 2*geom->Plane(plane)->Cell(0)->HalfL();
3423 // double cellW = 2*geom->Plane(plane)->Cell(0)->HalfW();
3424 // double cellD = 2*geom->Plane(plane)->Cell(0)->HalfD();
3425 // geo::View_t view = geom->Plane(plane)->View();
3426 /*
3427  if(view == geo::kX){
3428  std::cout<<"x view plane"<<std::endl;
3429  }else if(view == geo::kY){
3430  std::cout<<"y view plane"<<std::endl;
3431  }
3432  std::cout<<"cell L, cellW, cell D "<<cellL<<" "<<cellW<<" "<<cellD<<std::endl;
3433 */
3434  TVector3 toedge(xtoedge, ytoedge, ztoedge);
3435  return toedge;
3436  }
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 GetShwStop(rb::Prong shw)
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
float abs(float number)
Definition: d0nt_math.hpp:39
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
void geom(int which=0)
Definition: geom.C:163
TVector3 jmshower::RecoJMShower::GetTrkPlanePos ( unsigned int  plane,
rb::Prong  trk 
)

Definition at line 3133 of file RecoJMShower_module.cc.

References std::abs(), geo::PlaneGeo::Cell(), rb::Prong::Dir(), geom(), geo::CellGeo::GetCenter(), GetShwStop(), geo::GeometryBase::Plane(), and rb::Prong::Start().

Referenced by RecoShowers().

3133  {
3134  Double_t cellXYZ[3]={0,0,0};
3136  geom->Plane(plane)->Cell(0)->GetCenter(cellXYZ);
3137  TVector3 trkdir=trk.Dir();
3138  TVector3 trkstart=trk.Start();
3139  TVector3 trkstop=RecoJMShower::GetShwStop(trk);
3140  double trkplanex=99999,trkplaney=99999, trkplanez=0;
3141  if(std::abs(trkdir.z())>0.0001){
3142  trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3143  trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3144  }else{
3145  trkplanex = (trkstart.x()+trkstop.x())/2.0;
3146  trkplaney = (trkstart.y()+trkstop.y())/2.0;
3147  }
3148  trkplanez = cellXYZ[2];
3149  TVector3 trkplanepos(trkplanex, trkplaney, trkplanez);
3150  return trkplanepos;;
3151  }
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 GetShwStop(rb::Prong shw)
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
float abs(float number)
Definition: d0nt_math.hpp:39
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
void geom(int which=0)
Definition: geom.C:163
bool jmshower::RecoJMShower::IsFiducial ( rb::Prong  shower)

Definition at line 3439 of file RecoJMShower_module.cc.

References fShowerPClusterCol, GetTrkPlaneDistToEdge(), MECModelEnuComparisons::i, rb::Cluster::ID(), NDAPDHVSetting::plane, and Radius().

Referenced by RecoShowers().

3439  {
3440  int i = shower.ID();
3441  bool isfiducial = true;
3442  for(unsigned int plane=0;plane<fShowerPClusterCol[i].size();plane++){
3443  double planerad = RecoJMShower::Radius(fShowerPClusterCol[i][plane],shower);
3444  TVector3 disttoedge = RecoJMShower::GetTrkPlaneDistToEdge(fShowerPClusterCol[i][plane].Plane(), shower);
3445  if(disttoedge.x()>2*planerad&&disttoedge.y()>2*planerad&&disttoedge.z()>2*planerad)continue;
3446  isfiducial = false;
3447  }
3448  return isfiducial;
3449  }
std::vector< std::vector< jmshower::PCluster > > fShowerPClusterCol
double Radius(rb::Prong shower)
const int ID() const
Definition: Cluster.h:75
TVector3 GetTrkPlaneDistToEdge(unsigned int plane, art::Ptr< rb::Prong > trk)
bool jmshower::RecoJMShower::IsInvPhoton ( jmshower::JMShower  shower)

Definition at line 5255 of file RecoJMShower_module.cc.

References jmshower::JMShower::DedxInvLongLL(), jmshower::JMShower::DedxLongLL(), and jmshower::kGamma.

Referenced by RecoShowers().

5255  {
5256  if(shower.DedxLongLL(int(jmshower::kGamma))<shower.DedxInvLongLL(int(jmshower::kGamma))){
5257  return kTRUE;
5258  }else{
5259  return kFALSE;
5260  }
5261  }
double DedxInvLongLL(int type) const
Definition: JMShower.cxx:684
double DedxLongLL(int type) const
Definition: JMShower.cxx:672
bool jmshower::RecoJMShower::IsMuon ( jmshower::JMShower  shower)

Definition at line 5233 of file RecoJMShower_module.cc.

References jmshower::JMShower::DedxLongLL(), jmshower::JMShower::DedxTransLL(), makeTrainCVSamples::int, jmshower::kElectron, jmshower::kElectronCOH, jmshower::kElectronDIS, jmshower::kElectronQE, jmshower::kElectronRES, jmshower::kLast, jmshower::kMuon, and jmshower::kNull.

Referenced by RecoShowers().

5233  {
5234  double ellmax = -9999;
5235  double elllmax = -9999;
5236  double elltmax = -9999;
5237  for(int itype = int(jmshower::kNull)+1; itype != int(jmshower::kLast); ++itype){
5238  if(itype == int(jmshower::kElectron)||itype == int(jmshower::kElectronQE)||itype == int(jmshower::kElectronRES)||itype == int(jmshower::kElectronDIS)||itype == int(jmshower::kElectronCOH)){
5239  if(ellmax<shower.DedxLongLL(itype)+shower.DedxTransLL(itype)){
5240  ellmax = shower.DedxLongLL(itype)+shower.DedxTransLL(itype);
5241  elllmax = shower.DedxLongLL(itype);
5242  elltmax = shower.DedxTransLL(itype);
5243  }
5244  }
5245  }
5246  double emuLLL = elllmax - shower.DedxLongLL(int(jmshower::kMuon));
5247  double emuLLT = elltmax - shower.DedxTransLL(int(jmshower::kMuon));
5248  if(emuLLL>-0.2&&emuLLT>-99999){
5249  return kFALSE;
5250  }else{
5251  return kTRUE;
5252  }
5253  }
double DedxLongLL(int type) const
Definition: JMShower.cxx:672
double DedxTransLL(int type) const
Definition: JMShower.cxx:708
bool jmshower::RecoJMShower::LoadDedxHistograms ( int  detid)
private

Definition at line 488 of file RecoJMShower_module.cc.

References om::cout, allTimeWatchdog::endl, util::EnvExpansion(), fHtCellTransDedxE, fHtCellTransDedxEcoh, fHtCellTransDedxEdis, fHtCellTransDedxEqe, fHtCellTransDedxEres, fHtCellTransDedxEsg, fHtCellTransDedxG, fHtCellTransDedxMu, fHtCellTransDedxN, fHtCellTransDedxP, fHtCellTransDedxPi, fHtCellTransDedxPi0, fHtExpPlaneE, fHtExpPlaneEcoh, fHtExpPlaneEdis, fHtExpPlaneEqe, fHtExpPlaneEres, fHtExpPlaneEsg, fHtExpPlaneG, fHtExpPlaneHad, fHtExpPlaneMu, fHtExpPlaneN, fHtExpPlaneP, fHtExpPlanePi, fHtExpPlanePi0, fHtPlaneDedxE, fHtPlaneDedxEcoh, fHtPlaneDedxEdis, fHtPlaneDedxEqe, fHtPlaneDedxEres, fHtPlaneDedxEsg, fHtPlaneDedxG, fHtPlaneDedxMu, fHtPlaneDedxN, fHtPlaneDedxP, fHtPlaneDedxPi, fHtPlaneDedxPi0, fLibPath, GetEntries(), MECModelEnuComparisons::i, novadaq::cnv::kFARDET, novadaq::cnv::kNDOS, novadaq::cnv::kNEARDET, submit_hadd::l, pPIDHist, r(), and string.

Referenced by beginRun().

489  {
490 
492 
493  std::string fnameE(libpath);
494  std::string fnameG(libpath);
495  std::string fnameMu(libpath);
496  std::string fnamePi0(libpath);
497 // std::string fnameHad(libpath);
498  std::string fnameP(libpath);
499  std::string fnameN(libpath);
500  std::string fnamePi(libpath);
501  std::string fnameEqe(libpath);
502  std::string fnameEres(libpath);
503  std::string fnameEdis(libpath);
504  std::string fnameEcoh(libpath);
505  std::string fnameEsg(libpath);
506 
507  if(pPIDHist==1){
508  switch (detid) {
510 
511  fnameE += "e10kdedxhist_fd_nue_cc.root";
512  fnameG += "g10kdedxhist_fd_nue_cc.root";
513  fnameMu += "mu10kdedxhist_fd_nue_cc.root";
514  fnamePi0 += "pi010kdedxhist_fd_nue_cc.root";
515  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
516  fnameP += "protondedxhist_fd.root";
517  fnameN += "neutrondedxhist_fd.root";
518  fnamePi += "pidedxhist_fd.root";
519  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
520  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
521  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
522  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
523  fnameEsg += "singleehist_fd.root";
524 
525 /*
526  fnameE += "e10kdedxhist_fd_nue_cc.root";
527  fnameG += "g10kdedxhist_genie_fd_nhc.root";
528  fnameMu += "mu10kdedxhist_genie_fd_nhc.root";
529  fnamePi0 += "pi010kdedxhist_genie_fd_nhc.root";
530  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
531  fnameP += "protondedxhist_fd.root";
532  fnameN += "neutrondedxhist_fd.root";
533  fnamePi += "pidedxhist_fd.root";
534  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
535  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
536  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
537  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
538 */
539  break;
540  case novadaq::cnv::kNDOS:
541  fnameE += "e10kdedxhist_fd_nue_cc.root";
542  fnameG += "g10kdedxhist_fd_nue_cc.root";
543  fnameMu += "mu10kdedxhist_fd_nue_cc.root";
544  fnamePi0 += "pi010kdedxhist_fd_nue_cc.root";
545  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
546  fnameP += "protondedxhist_fd.root";
547  fnameN += "neutrondedxhist_fd.root";
548  fnamePi += "pidedxhist_fd.root";
549  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
550  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
551  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
552  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
553  fnameEsg += "singleehist_fd.root";
554  break;
556 
557  fnameE += "e10kdedxhist_fd_nue_cc.root";
558  fnameG += "g10kdedxhist_fd_nue_cc.root";
559  fnameMu += "mu10kdedxhist_fd_nue_cc.root";
560  fnamePi0 += "pi010kdedxhist_fd_nue_cc.root";
561  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
562  fnameP += "protondedxhist_fd.root";
563  fnameN += "neutrondedxhist_fd.root";
564  fnamePi += "pidedxhist_fd.root";
565  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
566  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
567  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
568  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
569  fnameEsg += "singleehist_fd.root";
570 
571  break;
572  default:
573  throw cet::exception("BAD DETID") << __FILE__ << ":" << __LINE__
574  << " Bad detector id = " << detid;
575  }
576  }else {
577  switch (detid) {
579 
580  fnameE += "e10kdedxhist_fd_nue_cc.root";
581  fnameG += "g10kdedxhist_fd_nue_cc.root";
582  fnameMu += "mu10kdedxhist_fd_nue_cc.root";
583  fnamePi0 += "pi010kdedxhist_fd_nue_cc.root";
584  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
585  fnameP += "protondedxhist_fd.root";
586  fnameN += "neutrondedxhist_fd.root";
587  fnamePi += "pidedxhist_fd.root";
588  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
589  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
590  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
591  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
592  fnameEsg += "singleehist_fd.root";
593 
594  break;
595  case novadaq::cnv::kNDOS:
596 
597  fnameE += "e10kdedxhist_fd_nue_cc.root";
598  fnameG += "g10kdedxhist_fd_nue_cc.root";
599  fnameMu += "mu10kdedxhist_fd_nue_cc.root";
600  fnamePi0 += "pi010kdedxhist_fd_nue_cc.root";
601  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
602  fnameP += "protondedxhist_fd.root";
603  fnameN += "neutrondedxhist_fd.root";
604  fnamePi += "pidedxhist_fd.root";
605  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
606  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
607  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
608  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
609  fnameEsg += "singleehist_fd.root";
610 
611  break;
613 
614  fnameE += "e10kdedxhist_fd_nue_cc.root";
615  fnameG += "g10kdedxhist_fd_nue_cc.root";
616  fnameMu += "mu10kdedxhist_fd_nue_cc.root";
617  fnamePi0 += "pi010kdedxhist_fd_nue_cc.root";
618  // fnameHad += "had10kdedxhist_fd_nue_ccqe.root";
619  fnameP += "protondedxhist_fd.root";
620  fnameN += "neutrondedxhist_fd.root";
621  fnamePi += "pidedxhist_fd.root";
622  fnameEqe += "e10kdedxhist_fd_nue_ccqe.root";
623  fnameEres += "e10kdedxhist_fd_nue_ccres.root";
624  fnameEdis += "e10kdedxhist_fd_nue_ccdis.root";
625  fnameEcoh += "e10kdedxhist_fd_nue_cccoh.root";
626  fnameEsg += "singleehist_fd.root";
627  break;
628  default:
629  throw cet::exception("BAD DETID") << __FILE__ << ":" << __LINE__
630  << " Bad detector id = " << detid;
631  }
632  }
633 
634 /* histograms stored in /nova/data/users/bianjm/share/development/RecoJMShower/histograms for the test versions
635 
636  std::string fFnameE;
637  std::string fFnameG;
638  std::string fFnameMu;
639  std::string fFnamePi0;
640  std::string fFnameHad;
641  std::string fFnameP;
642  std::string fFnameN;
643  std::string fFnamePi;
644  std::string fFnameEqe;
645  std::string fFnameEres;
646  std::string fFnameEdis;
647  std::string fFnameEcoh;
648 // constructor decides if initialized value is a path or an environment variable
649  cet::search_path sp("FW_SEARCH_PATH");
650  sp.find_file(fnameE, fFnameE);
651  sp.find_file(fnameG, fFnameG);
652  sp.find_file(fnameMu, fFnameMu);
653  sp.find_file(fnamePi0, fFnamePi0);
654  sp.find_file(fnameHad, fFnameHad);
655  sp.find_file(fnameP, fFnameP);
656  sp.find_file(fnameN, fFnameN);
657  sp.find_file(fnamePi, fFnamePi);
658  sp.find_file(fnameEqe, fFnameEqe);
659  sp.find_file(fnameEres, fFnameEres);
660  sp.find_file(fnameEdis, fFnameEdis);
661  sp.find_file(fnameEcoh, fFnameEcoh);
662 
663  struct stat sb;
664 
665  if (fFnameE.empty() || stat(fFnameE.c_str(), &sb)!=0){
666  // Failed to resolve the file name
667  mf::LogWarning("RecoJMShower") << "Histogram -- E " << fFnameE << " not found.";
668  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
669  }
670 
671  if (fFnameEqe.empty() || stat(fFnameEqe.c_str(), &sb)!=0){
672  // Failed to resolve the file name
673  mf::LogWarning("RecoJMShower") << "Histogram -- Eqe " << fFnameEqe << " not found.";
674  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
675  }
676  if (fFnameEres.empty() || stat(fFnameEres.c_str(), &sb)!=0){
677  // Failed to resolve the file name
678  mf::LogWarning("RecoJMShower") << "Histogram -- Eres " << fFnameEres << " not found.";
679  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
680  }
681  if (fFnameEdis.empty() || stat(fFnameEdis.c_str(), &sb)!=0){
682  // Failed to resolve the file name
683  mf::LogWarning("RecoJMShower") << "Histogram -- Edis " << fFnameEdis << " not found.";
684  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
685  }
686  if (fFnameEcoh.empty() || stat(fFnameEcoh.c_str(), &sb)!=0){
687  // Failed to resolve the file name
688  mf::LogWarning("RecoJMShower") << "Histogram -- Ecoh " << fFnameEcoh << " not found.";
689  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
690  }
691 
692  if (fFnameG.empty() || stat(fFnameG.c_str(), &sb)!=0){
693  // Failed to resolve the file name
694  mf::LogWarning("RecoJMShower") << "Histogram -- G " << fFnameG << " not found.";
695  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
696  }
697  if (fFnameMu.empty() || stat(fFnameMu.c_str(), &sb)!=0){
698  // Failed to resolve the file name
699  mf::LogWarning("RecoJMShower") << "Histogram -- Mu " << fFnameMu << " not found.";
700  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
701  }
702  if (fFnamePi0.empty() || stat(fFnamePi0.c_str(), &sb)!=0){
703  // Failed to resolve the file name
704  mf::LogWarning("RecoJMShower") << "Histogram -- Pi0 " << fFnamePi0 << " not found.";
705  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
706  }
707 
708  if (fFnameHad.empty() || stat(fFnameHad.c_str(), &sb)!=0){
709  // Failed to resolve the file name
710  mf::LogWarning("RecoJMShower") << "Histogram -- Had " << fFnameHad << " not found.";
711  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
712  }
713 
714  if (fFnameP.empty() || stat(fFnameP.c_str(), &sb)!=0){
715  // Failed to resolve the file name
716  mf::LogWarning("RecoJMShower") << "Histogram -- P " << fFnameP << " not found.";
717  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
718  }
719 
720  if (fFnameN.empty() || stat(fFnameN.c_str(), &sb)!=0){
721  // Failed to resolve the file name
722  mf::LogWarning("RecoJMShower") << "Histogram -- N " << fFnameN << " not found.";
723  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
724  }
725 
726  if (fFnamePi.empty() || stat(fFnamePi.c_str(), &sb)!=0){
727  // Failed to resolve the file name
728  mf::LogWarning("RecoJMShower") << "Histogram -- Pi " << fFnamePi << " not found.";
729  throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
730  }
731 
732 */
733  TFile hfileE(fnameE.c_str());
734  TFile hfileG(fnameG.c_str());
735  TFile hfileMu(fnameMu.c_str());
736  TFile hfilePi0(fnamePi0.c_str());
737 // TFile hfileHad(fnameHad.c_str());
738  TFile hfileP(fnameP.c_str());
739  TFile hfileN(fnameN.c_str());
740  TFile hfilePi(fnamePi.c_str());
741  TFile hfileEqe(fnameEqe.c_str());
742  TFile hfileEres(fnameEres.c_str());
743  TFile hfileEdis(fnameEdis.c_str());
744  TFile hfileEcoh(fnameEcoh.c_str());
745  TFile hfileEsg(fnameEsg.c_str());
746  std::cout<<"PID Dedx File "<<fnameE.c_str()<<" is loaded."<<std::endl;
747  std::cout<<"PID Dedx File "<<fnameG.c_str()<<" is loaded."<<std::endl;
748  std::cout<<"PID Dedx File "<<fnameMu.c_str()<<" is loaded."<<std::endl;
749  std::cout<<"PID Dedx File "<<fnamePi0.c_str()<<" is loaded."<<std::endl;
750 // std::cout<<"PID Dedx File "<<fnameHad.c_str()<<" is loaded."<<std::endl;
751  std::cout<<"PID Dedx File "<<fnameP.c_str()<<" is loaded."<<std::endl;
752  std::cout<<"PID Dedx File "<<fnameN.c_str()<<" is loaded."<<std::endl;
753  std::cout<<"PID Dedx File "<<fnamePi.c_str()<<" is loaded."<<std::endl;
754  std::cout<<"PID Dedx File "<<fnameEqe.c_str()<<" is loaded."<<std::endl;
755  std::cout<<"PID Dedx File "<<fnameEres.c_str()<<" is loaded."<<std::endl;
756  std::cout<<"PID Dedx File "<<fnameEdis.c_str()<<" is loaded."<<std::endl;
757  std::cout<<"PID Dedx File "<<fnameEcoh.c_str()<<" is loaded."<<std::endl;
758  std::cout<<"PID Dedx File "<<fnameEsg.c_str()<<" is loaded."<<std::endl;
759 
760  //The detector is divided into 4 regions in XY, these regions extend in z.
761  //Region 1 is the top left quadrant in XY when looking at the front of the detector
762  //The training histograms are then divided into 11 energy regions from 0 to 5 GeV.
763  //In the training each shower is put into a bin based on its reco energy
764  //The first and last energy bin are 0-0.25 and 4.75-5 GeV. All bins inbetween are 0.5 GeV wide.
765  //For each region and each energy there is a de/dx histogram for each plane from 0 to 200
766  //the histograms being loaded below are for longitudinal dedx.
767  //These training histograms are also broken up by particle type (electron, muon, gamma, proton, pi0, pion, neutron)
768 
769  //NOTE: Jianming is exploring using the shower length as a training variable.
770  //In order to do this the variables below are initialized and then store the
771  //shower length at each energy in each region. This is done by finding the
772  //last plane in which there was energy deposited.
773  //These variables are not currently being used in the production version of the algorithm
774  for(int r=0; r<4; r++){
775  for(int i=0; i<11; i++){
776  fHtExpPlaneE[r][i]=0;
777  fHtExpPlaneG[r][i]=0;
778  fHtExpPlaneMu[r][i]=0;
779  fHtExpPlanePi0[r][i]=0;
780  fHtExpPlaneHad[r][i]=0;
781  fHtExpPlaneP[r][i]=0;
782  fHtExpPlaneN[r][i]=0;
783  fHtExpPlanePi[r][i]=0;
784  fHtExpPlaneEqe[r][i]=0;
785  fHtExpPlaneEres[r][i]=0;
786  fHtExpPlaneEdis[r][i]=0;
787  fHtExpPlaneEcoh[r][i]=0;
788  }
789  }
790 
791  for(int r=0; r<4; r++){
792  for(int i=0; i<41; i++){
793  fHtExpPlaneEsg[r][i]=0;
794  }
795  }
796 
797 
798 
799  for(int r=0; r<4; r++){
800  for(int i=0; i<11; i++){
801  for(int l=0; l<200; l++){
802  char cht[200];
803  sprintf(cht,"ht_%d_%d_%d",r,i,l);
804  fHtPlaneDedxE[r][i][l] = (TH1D*)hfileE.Get(cht);
805  fHtPlaneDedxG[r][i][l] = (TH1D*)hfileG.Get(cht);
806  fHtPlaneDedxMu[r][i][l] = (TH1D*)hfileMu.Get(cht);
807  fHtPlaneDedxPi0[r][i][l] = (TH1D*)hfilePi0.Get(cht);
808  fHtPlaneDedxP[r][i][l] = (TH1D*)hfileP.Get(cht);
809  fHtPlaneDedxN[r][i][l] = (TH1D*)hfileN.Get(cht);
810  fHtPlaneDedxPi[r][i][l] = (TH1D*)hfilePi.Get(cht);
811  fHtPlaneDedxEqe[r][i][l] = (TH1D*)hfileEqe.Get(cht);
812  fHtPlaneDedxEres[r][i][l] = (TH1D*)hfileEres.Get(cht);
813  fHtPlaneDedxEdis[r][i][l] = (TH1D*)hfileEdis.Get(cht);
814  fHtPlaneDedxEcoh[r][i][l] = (TH1D*)hfileEcoh.Get(cht);
815 
816 
817  if(fHtPlaneDedxE[r][i][l]->GetEntries()>0&&fHtExpPlaneE[r][i]==0){
818  if(fHtPlaneDedxE[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxE[r][i][l]->GetEntries()>0.997){
819  fHtExpPlaneE[r][i]=l;
820  }
821  }
822 
823  if(fHtPlaneDedxG[r][i][l]->GetEntries()>0&&fHtExpPlaneG[r][i]==0){
824  if(fHtPlaneDedxG[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxG[r][i][l]->GetEntries()>0.997){
825  fHtExpPlaneG[r][i]=l;
826  }
827  }
828 
829  if(fHtPlaneDedxMu[r][i][l]->GetEntries()>0&&fHtExpPlaneMu[r][i]==0){
830  if(fHtPlaneDedxMu[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxMu[r][i][l]->GetEntries()>0.997){
831  fHtExpPlaneMu[r][i]=l;
832  }
833  }
834 
835  if(fHtPlaneDedxPi0[r][i][l]->GetEntries()>0&&fHtExpPlanePi0[r][i]==0){
836  if(fHtPlaneDedxPi0[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxPi0[r][i][l]->GetEntries()>0.997){
837  fHtExpPlanePi0[r][i]=l;
838  }
839  }
840 
841  if(fHtPlaneDedxP[r][i][l]->GetEntries()>0&&fHtExpPlaneP[r][i]==0){
842  if(fHtPlaneDedxP[r][i][l]->GetBinContent(1)/(double)fHtPlaneDedxP[