LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
LinearEnergyAlg.cxx
Go to the documentation of this file.
1 
10 // LArSoft libraries
12 #include "lardata/Utilities/ForEachAssociatedGroup.h" // util::associated_groups()
14 
15 // Framework libraries
17 
18 // C/C++ standard libraries
19 #include <stdexcept> // std::runtime_error()
20 #include <cmath>
21 
22 
23 namespace {
24 
32  auto hitsAssociatedWith(
34  art::Assns<recob::Cluster, recob::Hit> const& hitsPerCluster
35  )
36  {
37  //
38  // Reminder: the cluster-hit association is a list of links:
39  //
40  // cluster #0 <--> hit #7
41  // cluster #0 <--> hit #8
42  // cluster #0 <--> hit #10
43  // cluster #1 <--> hit #4
44  // cluster #1 <--> hit #5
45  // cluster #1 <--> hit #9
46  // cluster #2 <--> hit #11
47  // cluster #2 <--> hit #12
48  // ...
49  //
50  // associated_groups() presents this associations grouped by cluster:
51  //
52  // cluster #0 <--> hits { #7, #8, #10 }
53  // cluster #1 <--> hits { #4, #5, #9 }
54  // cluster #2 <--> hits { #11, #12 }
55  // ...
56  //
57  // This function loops through this latter list, looking for the cluster
58  // number that is stored in the `cluster` argument.
59  //
60  size_t iAssociatedCluster = 0;
61  for (auto hits: util::associated_groups(hitsPerCluster)) {
62 
63  // if the iAssociatedCluster is the index we are looking for, we are done:
64  if (iAssociatedCluster == cluster.key()) return hits;
65 
66  ++iAssociatedCluster; // else, try the next one!
67 
68  } // for all clusters
69  throw std::runtime_error
70  ("No hit associated with cluster " + std::to_string(cluster.key()) + " found!");
71  } // calo::LinearEnergyAlg::hitsAssociatedWith()
72 
73 } // local namespace
74 
75 
76 //------------------------------------------------------------------------------
77 const std::string calo::LinearEnergyAlg::ModelName::ModBox = "ModBox";
78 const std::string calo::LinearEnergyAlg::ModelName::Birks = "Birks";
79 const std::string calo::LinearEnergyAlg::ModelName::Constant = "Constant";
80 
81 
83  : fUseArea( config.UseArea() )
84 // , fRecombFactor( 1. )
85  , fElectronLifetime( 1e10 ) // needs to be read from service
86  , fDeconNorm( 200 )
87 {
88  auto const& recombParams = config.Recombination();
89  if ( recombParams.modelIsModBox()) {
91  recombModBoxParams.A = recombParams.A();
92  recombModBoxParams.B = recombParams.B();
93  }
94  else if ( recombParams.modelIsBirks()) {
96  recombBirksParams.A = recombParams.A3t();
97  recombBirksParams.k = recombParams.k3t();
98  }
99  else if ( recombParams.modelIsConstant()) {
101  recombConstParams.factor = recombParams.factor();
102  }
103  else {
104  throw std::runtime_error
105  ( "Unsupported recombination mode: '" + recombParams.Model() + "'" );
106  }
107 
108  if (!fUseArea) {
109  throw std::runtime_error("Energy correction based on hit peak amplitude not implemented yet.");
110  }
111 
112 }
113 
114 
116 
118 
119 }
120 
121 
122 
124 {
125 
126  double const t = detc->TPCTick2TrigTime( hit.PeakTime() );
127  double const LifetimeCorr = std::exp( t / fElectronLifetime );
128 
129  // hit charge (ADC) -> Coulomb -> Number of electrons -> eV
130  // TODO: implement the non-UseArea option
131  double dE = hit.Integral() * kWion * fDeconNorm;
132 
133  // dE * lifetime correction
134  dE *= LifetimeCorr;
135 
136  // dE / recombination factor
137  double const dEdx = 2.3;
138  double const RecombCorr = RecombinationCorrection( dEdx ) * kWion / dEdx;
139  dE /= RecombCorr;
140 
141  return dE;
142 
143 } // calo::LinearEnergyAlg::CalculateHitEnergy()
144 
145 
146 
148  std::vector<art::Ptr<recob::Cluster>> const& clusters,
149  art::Assns<recob::Cluster, recob::Hit> const& hitsPerCluster
150  ) const
151 { // input clusters and hits, shower direction?
152 
153  if (clusters.empty()) return {}; // no clusters!!
154 
155  // prepare the energies vector with one entry per plane
156  // (we get the total number of planes of the TPC the cluster is in from geometry)
157  // and initialize them to a ridiculously negative number to start with
158 
159  std::vector<double> clusterEnergies;
160  geo::TPCID refTPC = clusters[0]->Plane();
161  clusterEnergies.resize
162  (geom->TPC(refTPC).Nplanes(), std::numeric_limits<double>::lowest());
163 
164  if ( clusters.size() > clusterEnergies.size() ) {
165  mf::LogError("LinearEnergyAlg") << clusters.size() << " clusters for "
166  << clusterEnergies.size() << " wire planes!";
167  }
168 
169  for (art::Ptr<recob::Cluster> const& cluster: clusters) {
170 
171  auto const plane = cluster->Plane();
172  if (plane != refTPC) {
173  throw std::runtime_error(
174  "Cluster ID=" + std::to_string(cluster->ID())
175  + " is expected on TPC " + std::string(refTPC)
176  + " but is found on plane " + std::string(plane)
177  );
178  }
179 
180  // hitsAssociatedWith() searches for the right association links
181  // to the cluster we are processing
182  double const E = CalculateClusterEnergy
183  (*cluster, hitsAssociatedWith(cluster, hitsPerCluster));
184 
185  auto const planeNo = plane.Plane;
186  if (clusterEnergies[planeNo] >= 0.) {
187  mf::LogWarning("LinearEnergyAlg")
188  << "Warning! two or more clusters share plane "
189  << plane << "! (the last with energy " << E
190  << ", the previous " << clusterEnergies[planeNo] << " GeV)";
191  }
192  clusterEnergies[planeNo] = E;
193 
194  } // for all clusters
195 
196  return clusterEnergies;
197 
198 } // calo::LinearEnergyAlg::CalculateEnergy()
199 
201 
202  switch ( fRecombModel ) {
203  case kModBox:
204  return this->ModBoxInverse( dEdx );
205  case kBirks:
206  return this->BirksInverse( dEdx );
207  case kConstant:
208  return recombConstParams.factor;
209  default:
210  throw std::logic_error
211  ("Unexpected! recombination model in RecombinationCorrection()");
212  }
213 }
214 
215 // Modified Box model correction, should be moved to somewhere else in the future
216 double calo::LinearEnergyAlg::ModBoxInverse( double dEdx ) const {
217  // Modified Box model correction has better behavior than the Birks
218  // correction at high values of dQ/dx.
219  double rho = detp->Density(); // LAr density in g/cm^3
220  double Efield = detp->Efield(); // Electric Field in the drift region in KV/cm
221  double Beta = recombModBoxParams.B / (rho * Efield);
222  double Alpha = recombModBoxParams.A;
223 
224  double dQdx = std::log ( Alpha + Beta * dEdx ) / ( Beta * kWion );
225 
226  return dQdx;
227 }
228 
229 double calo::LinearEnergyAlg::BirksInverse( double dEdx ) const {
230  // Correction for charge quenching using parameterization from
231  // S.Amoruso et al., NIM A 523 (2004) 275
232 
233  double A3t = recombBirksParams.A;
234  double K3t = recombBirksParams.k; // in KV/cm*(g/cm^2)/MeV
235  double rho = detp->Density(); // LAr density in g/cm^3
236  double Efield = detp->Efield(); // Electric Field in the drift region in KV/cm
237  K3t /= rho; // KV/MeV
238 
239  double dQdx = ( A3t/kWion ) / ( K3t / Efield * dEdx + 1);
240 
241  return dQdx;
242 }
key_type key() const
Definition: Ptr.h:356
static const std::string Constant
geo::GeometryCore const * geom
Pointer to the geometry to be used.
Float_t E
Definition: plot.C:23
Algorithm(s) calculating energy.
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:145
std::vector< double > CalculateEnergy(std::vector< art::Ptr< recob::Cluster >> const &clusters, art::Assns< recob::Cluster, recob::Hit > const &hitsPerCluster) const
Calculates the energy of the shower.
double RecombinationCorrection(double dEdx) const
TODO: make it more flexible.
fhicl::Table< RecombinationConfig > Recombination
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:225
geo::PlaneID Plane() const
Returns the plane ID this cluster lies on.
Definition: Cluster.h:744
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Cluster finding and building.
double BirksInverse(double dEdx) const
Algorithm configuration.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
auto associated_groups(A const &assns)
Helper functions to access associations in order.
static const std::string Birks
void hits()
Definition: readHits.C:15
double CalculateClusterEnergy(recob::Cluster const &cluster, BeginHitIter beginHit, EndHitIter endHit) const
Calculates the energy of a single cluster.
double ModBoxInverse(double dEdx) const
detinfo::DetectorClocks const * detc
Pointer to the detector clock.
static constexpr double kWion
ionization potenial in LAr, 23.6 eV = 1e, Wion in GeV/e
The data type to uniquely identify a TPC.
Definition: geo_types.h:195
ModBoxParameters recombModBoxParams
Parameters for recombination box model; filled only when this model is selected.
virtual double Density(double temperature) const =0
Returns argon density at a given temperature.
Definition of data types for geometry description.
detinfo::DetectorProperties const * detp
Pointer to the detector property.
Detector simulation of raw signals on wires.
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:219
ConstantRecombParameters recombConstParams
Parameters for constant recombination factor; filled only when this model is selected.
std::string to_string(Flag_t< Storage > const flag)
Convert a flag into a stream (shows its index).
Definition: BitMask.h:187
ID_t ID() const
Identifier of this cluster.
Definition: Cluster.h:738
Helper functions to access associations in order.
LinearEnergyAlg(Config const &config)
Constructor with configuration validation.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
BirksParameters recombBirksParams
Parameters for recombination Birks model; filled only when this model is selected.
static const std::string ModBox
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
virtual double TPCTick2TrigTime(double tick) const =0
Converts a TPC time (in ticks) into a trigger time [µs].
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:49
virtual double Efield(unsigned int planegap=0) const =0
Returns the nominal electric field in the specified volume.
virtual double ElectronLifetime() const =0
Returns the attenuation constant for ionization electrons.
double CalculateHitEnergy(recob::Hit const &hit) const
Returns the corrected energy from the hit.