ProtoDUNEShowerUtils.cxx
Go to the documentation of this file.
2 
6 #include "canvas/Persistency/Common/FindManyP.h"
7 
9 
10 }
11 
13 
14 }
15 
16 // Get the hits from a given reco shower
17 const std::vector<const recob::Hit*> protoana::ProtoDUNEShowerUtils::GetRecoShowerHits(const recob::Shower &shower, art::Event const &evt, const std::string showerModule) const{
18 
19  auto recoShowers = evt.getValidHandle<std::vector<recob::Shower> >(showerModule);
20  art::FindManyP<recob::Hit> findHits(recoShowers,evt,showerModule);
21 
22  // Shower.ID is sometimes at a default value - make sure we get the correct one
23  int actualIndex = GetShowerIndex(shower,evt,showerModule);
24 
25  std::vector<art::Ptr<recob::Hit>> inputHits = findHits.at(actualIndex);
26  std::vector<const recob::Hit*> showerHits;
27 
28  for(const art::Ptr<recob::Hit> hit : inputHits){
29 
30  showerHits.push_back(hit.get());
31 
32  }
33 
34  return showerHits;
35 
36 }
37 
38 
39 // Get the hits from a given reco shower
40 const std::vector<art::Ptr<recob::Hit>> protoana::ProtoDUNEShowerUtils::GetRecoShowerArtHits(
41  const recob::Shower &shower, art::Event const &evt, const std::string showerModule) const {
42 
43  auto recoShowers = evt.getValidHandle<std::vector<recob::Shower> >(showerModule);
44  art::FindManyP<recob::Hit> findHits(recoShowers,evt,showerModule);
45 
46  // Shower.ID is sometimes at a default value - make sure we get the correct one
47  int actualIndex = GetShowerIndex(shower,evt,showerModule);
48 
49  std::vector<art::Ptr<recob::Hit>> showerHits = findHits.at(actualIndex);
50  return showerHits;
51 
52 }
53 
54 // Get the hits from a given reco shower
56 
57  return GetRecoShowerHits(shower,evt,showerModule).size();
58 
59 }
60 
61 // Get the PCAxis object from the reco shower
62 std::vector<const recob::PCAxis*> protoana::ProtoDUNEShowerUtils::GetRecoShowerPCAxis(const recob::Shower &shower, art::Event const &evt, const std::string showerModule) const{
63 
64  auto recoShowers = evt.getValidHandle<std::vector<recob::Shower> >(showerModule);
65  art::FindManyP<recob::PCAxis> findPCA(recoShowers,evt,showerModule);
66 
67  // Shower.ID is sometimes at a default value - make sure we get the correct one
68  int actualIndex = GetShowerIndex(shower,evt,showerModule);
69 
70  std::vector<const recob::PCAxis*> pcaVec;
71  for(auto const pca : findPCA.at(actualIndex)){
72  pcaVec.push_back(pca.get());
73  }
74 
75  return pcaVec;
76 }
77 
79  detinfo::DetectorPropertiesData const& detProp,
80  const std::vector<const recob::Hit*> &hits, calo::CalorimetryAlg caloAlg)
81 {
82  double kGeVtoElectrons { 4.237e7 }; // obtained from utils class.. Copied for now, should use class (although this is a physical constant, so hopefully doesn't change).
83  double recombination { 1/0.63 };
84 
85  std::vector<double> showerEnergy = {0,0,0};
86 
87  // Find the total charge on each plane
88  for ( size_t h{0} ; h < hits.size() ; h++ ) {
89  const recob::Hit* hit = hits[h];
90  const int plane = hit->WireID().Plane;
91  showerEnergy[ plane ] += ( caloAlg.ElectronsFromADCArea( hit->Integral(), plane) * caloAlg.LifetimeCorrection(clockData, detProp, hit->PeakTime()) ) / kGeVtoElectrons;
92  }
93 
94  showerEnergy[0] *= recombination;
95  showerEnergy[1] *= recombination;
96  showerEnergy[2] *= recombination;
97 
98  // caloAlg.ElectronsFromADCArea( hit->Integral(), plane) -> Does hit->Integral()/AreaConstants(plane)
99  // AreaConstants(plane) is defined in calorimetry_pdune.fcl. Although that fcl file has a typo.
100  // These probably need tuning for protodune data.
101 
102  return showerEnergy;
103 }
104 
105 // If the shower.ID() isn't filled we must find the actual shower index ourselves
107 
108  if(shower.ID() != -999) return shower.ID();
109 
110  auto recoShowers = evt.getValidHandle<std::vector<recob::Shower> >(showerModule);
111 
112  // Iterate through all showers to find the matching one to our shower
113  int actualIndex = shower.ID();
114  if(shower.ID() < 0){
115  for(unsigned int s = 0; s < recoShowers->size(); ++s){
116  const recob::Shower thisShower = (*recoShowers)[s];
117  // Can't compare actual objects so look at a property
118  if(fabs(thisShower.Length() - shower.Length()) < 1e-5){
119  actualIndex = s;
120  continue;
121  }
122  }
123  }
124 
125  return actualIndex;
126 
127 }
128 
129 
130 std::vector<anab::Calorimetry> protoana::ProtoDUNEShowerUtils::GetRecoShowerCalorimetry(const recob::Shower &shower, art::Event const &evt, const std::string showerModule, const std::string caloModule) const{
131 
132  auto recoShowers = evt.getValidHandle<std::vector<recob::Shower> >(showerModule);
133  std::vector<anab::Calorimetry> caloInfo;
134 
135  try{
136  const art::FindManyP<anab::Calorimetry> findCalorimetry(recoShowers,evt,caloModule);
137  int actualIndex = GetShowerIndex(shower,evt,showerModule);
138  std::vector<art::Ptr<anab::Calorimetry>> theseCalos = findCalorimetry.at(actualIndex);
139 
140  for( auto calo : theseCalos){
141  caloInfo.push_back(*calo);
142  }
143  }
144  catch(...){
145  std::cerr << "No shower calorimetry object found... returning empty vector" << std::endl;
146  }
147 
148  return caloInfo;
149 
150 
151 }
152 
153 // Get the space points associated to the Shower
154 /*
155 const std::vector<const recob::SpacePoint*>
156  protoana::ProtoDUNEShowerUtils::GetShowerSpacePoints(
157  const recob::Shower & shower, art::Event const & evt,
158  const std::string showerLabel) const {
159 
160  // Get the particles and their associations
161  auto particles =
162  evt.getValidHandle<std::vector<recob::PFParticle>>(showerLabel);
163  const art::FindManyP<recob::SpacePoint>
164  findSpacePoints(particles, evt, showerLabel);
165  const std::vector<art::Ptr<recob::SpacePoint>> pfpSpacePoints = findSpacePoints.at(particle.Self());
166 
167  // We don't want the art::Ptr so we need to get rid of it
168  std::vector<const recob::SpacePoint*> sp;
169  for(auto pointer : pfpSpacePoints){
170  sp.push_back(pointer.get());
171  }
172 
173  return sp;
174 }
175 */
const std::vector< const recob::Hit * > GetRecoShowerHits(const recob::Shower &shower, art::Event const &evt, const std::string showerModule) const
Get the hits from a given reco shower.
std::vector< double > EstimateEnergyFromHitCharge(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const std::vector< const recob::Hit * > &hits, calo::CalorimetryAlg caloAlg)
double Length() const
Definition: Shower.h:201
std::string string
Definition: nybbler.cc:12
geo::WireID WireID() const
Definition: Hit.h:233
std::vector< anab::Calorimetry > GetRecoShowerCalorimetry(const recob::Shower &shower, art::Event const &evt, const std::string showerModule, const std::string caloModule) const
Get shower calo info.
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
uint size() const
Definition: qcstring.h:201
unsigned int GetNumberRecoShowerHits(const recob::Shower &shower, art::Event const &evt, const std::string showerModule) const
Get the number of hits from a given reco shower.
double ElectronsFromADCArea(double area, unsigned short plane) const
const std::vector< art::Ptr< recob::Hit > > GetRecoShowerArtHits(const recob::Shower &shower, art::Event const &evt, const std::string showerModule) const
const double e
int GetShowerIndex(const recob::Shower &shower, art::Event const &evt, const std::string showerModule) const
If the shower.ID() isn&#39;t filled we must find the actual shower index ourselves.
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
double recombination
Definition: doAna.cpp:19
std::vector< const recob::PCAxis * > GetRecoShowerPCAxis(const recob::Shower &shower, art::Event const &evt, const std::string showerModule) const
Get the associated PCAxis object (from a principal component analysis)
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
Detector simulation of raw signals on wires.
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
Declaration of signal hit object.
Contains all timing reference information for the detector.
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
TCEvent evt
Definition: DataStructs.cxx:7
double LifetimeCorrection(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, double time, double T0=0) const
int ID() const
Definition: Shower.h:187
static QCString * s
Definition: config.cpp:1042
calorimetry
QTextStream & endl(QTextStream &s)