EnergyLossVsDistance.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file EnergyLossVsDistance.cxx
3 // \brief Integrate energy losses to calculate deposit near a track end
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
8 
9 #include <algorithm>
10 
11 #include "TF1.h"
12 #include "TGraph.h"
13 
14 namespace calib
15 {
16 
17  //......................................................................
19  double localDensity,
20  double bulkDensity)
21  {
22  // Start at slightly nonzero energy
23  double T = 2;
24  // up to 100m in 1mm increments
25  const double dx = .1;
26  for(double x = 0; x < 1e6; x += dx){
27  const double dEdx = bb->dEdx(T);
28 
29  fX_pts.push_back(x_pt(x, localDensity*dEdx));
30 
31  T += dEdx*dx*bulkDensity;
32  }
33  }
34 
36  {
37  return pt.x < x;
38  }
39  //......................................................................
41  {
42  typedef std::vector<x_pt>::const_iterator it_t;
43 
44  // Find the point after x
45  it_t itNext = std::lower_bound(fX_pts.begin(), fX_pts.end(), x, ltXPt);
46  if(itNext == fX_pts.end()) --itNext;
47  if(itNext == fX_pts.begin()) ++itNext;
48  // And the point before
49  it_t itPrev = itNext; --itPrev;
50 
51  // Interpolate
52  const double x0 = itPrev->x; const double y0 = itPrev->dEdx;
53  const double x1 = itNext->x; const double y1 = itNext->dEdx;
54  return ((x1-x)*y0+(x-x0)*y1)/(x1-x0);
55  }
56 
57  //......................................................................
59  {
60  TGraph* ret = new TGraph;
61  // 1mm to 100m
62  for(double x = .1; x < 1e4; x *= 1.05){
63  ret->SetPoint(ret->GetN(), x, GetEnergyLoss(x));
64  }
65 
66  return ret;
67  }
68 
69  //......................................................................
70  double EnergyLossVsDistance::Eval(double* xs, double* pars)
71  {
72  const double x = xs[0];
73  const double scale = pars[0];
74 
75  return scale*GetEnergyLoss(x);
76  }
77 
78  //......................................................................
80  {
81  TF1* ret = new TF1("el", this, &EnergyLossVsDistance::Eval,
82  0, 500, 1);
83  ret->SetParName(0, "energy_scale");
84  ret->SetParameter(0, 1);
85  return ret;
86  }
87 
88 } // end namespace
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
double GetEnergyLoss(double dist) const
virtual double dEdx(double T) const =0
bool ltXPt(EnergyLossVsDistance::x_pt pt, double x)
Double_t scale
Definition: plot.C:25
CDPStorage service.
double dx[NP][NC]
EnergyLossVsDistance(IBetheBloch *bb, double localDensity, double bulkDensity)
std::string pars("Th23Dmsq32")
double T
Definition: Xdiff_gwt.C:5
double Eval(double *xs, double *pars)