CbmRoot
CbmStsPhysics.cxx
Go to the documentation of this file.
1 
6 #include "CbmStsPhysics.h"
7 
8 #include "CbmStsDefs.h" // for CbmSts, kProtonMass, kSiCharge, kSiDe...
9 
10 #include <FairLogger.h> // for Logger, LOG
11 
12 #include <TDatabasePDG.h> // for TDatabasePDG
13 #include <TMath.h> // for Log, Power
14 #include <TParticlePDG.h> // for TParticlePDG
15 #include <TRandom.h> // for TRandom, gRandom
16 #include <TString.h> // for TString, operator<<, operator+
17 #include <TSystem.h> // for TSystem, gSystem
18 
19 #include <fstream> // for ifstream, basic_istream, right
20 #include <iomanip> // for setw, __iom_t6
21 #include <math.h> // for log, sqrt
22 #include <utility> // for pair
23 
24 #include "CbmStsPhysics.h"
25 
26 using std::ifstream;
27 using std::map;
28 using std::right;
29 using std::setw;
30 using namespace CbmSts;
31 
33 
34  // ----- Initialisation of static variables ----------------------------
36 // -------------------------------------------------------------------------
37 
38 
39 // ----- Constructor ---------------------------------------------------
41  // --- Read the energy loss data tables
42  LOG(info) << "Instantiating STS Physics... ";
43  ReadDataTablesStoppingPower();
44  ReadDataTablesLandauWidth();
45 
46  // --- Initialise the constants for the Urban model
47  SetUrbanParameters(CbmSts::kSiCharge);
48 }
49 // -------------------------------------------------------------------------
50 
51 
52 // ----- Destructor ----------------------------------------------------
54 // -------------------------------------------------------------------------
55 
56 
57 // ----- Diffusion width -----------------------------------------------
58 Double_t CbmStsPhysics::DiffusionWidth(Double_t z,
59  Double_t d,
60  Double_t vBias,
61  Double_t vFd,
62  Double_t temperature,
63  Int_t chargeType) {
64 
65  // --- Check parameters. A tolerance of 0.1 micrometer on the sensor borders
66  // --- is used to avoid crashes due to rounding errors.
67  if (z < 0. && z > -0.00001) z = 0.;
68  if (z > d && z < d + 0.00001) z = d;
69  if (z < 0. || z > d) {
70  LOG(error) << "StsPhysics: z coordinate " << z
71  << " not inside sensor (d = " << d << ")";
72  return -1.;
73  }
74  if (temperature < 0.) {
75  LOG(error) << "StsPhysics: illegal temperature value " << temperature;
76  return -1.;
77  }
78 
79  // --- Diffusion constant over mobility [J/C]
80  // --- The numerical factor is k_B/e in units of J/(KC).
81  Double_t diffConst = 8.61733e-5 * temperature;
82 
83  // --- Drift time times mobility [cm**2 * C / J]
84  // For the formula, see the STS digitiser note.
85  Double_t tau = 0.;
86  if (chargeType == 0) { // electrons, drift to n (front) side
87  tau = 0.5 * d * d / vFd
88  * log((vBias + (1. - 2. * z / d) * vFd) / (vBias - vFd));
89  } else if (chargeType == 1) { // holes, drift to the p (back) side
90  tau = -0.5 * d * d / vFd * log(1. - 2. * vFd * z / d / (vBias + vFd));
91  } else {
92  LOG(error) << "StsPhysics: Illegal charge type " << chargeType;
93  return -1.;
94  }
95 
96  return sqrt(2. * diffConst * tau);
97 }
98 // -------------------------------------------------------------------------
99 
100 
101 // ----- Electric field ------------------------------------------------
102 Double_t CbmStsPhysics::ElectricField(Double_t vBias,
103  Double_t vFd,
104  Double_t dZ,
105  Double_t z) {
106  return (vBias + vFd * (2. * z / dZ - 1.)) / dZ;
107 }
108 // -------------------------------------------------------------------------
109 
110 
111 // ----- Energy loss from fluctuation model ----------------------------
112 Double_t CbmStsPhysics::EnergyLoss(Double_t dz,
113  Double_t mass,
114  Double_t eKin,
115  Double_t dedx) const {
116 
117  // Gamma and beta
118  Double_t gamma = (eKin + mass) / mass;
119  Double_t beta2 = 1. - 1. / (gamma * gamma);
120 
121  // Auxiliary
122  Double_t xAux = 2. * mass * beta2 * gamma * gamma;
123 
124  // Mean energy losses (PHYS333 2.4 eqs. (2) and (3))
125  Double_t sigma1 = dedx * fUrbanF1 / fUrbanE1
126  * (TMath::Log(xAux / fUrbanE1) - beta2)
127  / (TMath::Log(xAux / fUrbanI) - beta2) * (1. - fUrbanR);
128  Double_t sigma2 = dedx * fUrbanF2 / fUrbanE2
129  * (TMath::Log(xAux / fUrbanE2) - beta2)
130  / (TMath::Log(xAux / fUrbanI) - beta2) * (1. - fUrbanR);
131  Double_t sigma3 = dedx * fUrbanEmax * fUrbanR
132  / (fUrbanI * (fUrbanEmax + fUrbanI))
133  / TMath::Log((fUrbanEmax + fUrbanI) / fUrbanI);
134 
135  // Sample number of processes Poissonian energy loss distribution
136  // (PHYS333 2.4 eq. (6))
137  Int_t n1 = gRandom->Poisson(sigma1 * dz);
138  Int_t n2 = gRandom->Poisson(sigma2 * dz);
139  Int_t n3 = gRandom->Poisson(sigma3 * dz);
140 
141  // Ion energy loss (PHYS333 2.4 eq. (12))
142  Double_t eLossIon = 0.;
143  for (Int_t j = 1; j <= n3; j++) {
144  Double_t uni = gRandom->Uniform(1.);
145  eLossIon += fUrbanI / (1. - uni * fUrbanEmax / (fUrbanEmax + fUrbanI));
146  }
147 
148  // Total energy loss
149  return (n1 * fUrbanE1 + n2 * fUrbanE2 + eLossIon);
150 }
151 // -------------------------------------------------------------------------
152 
153 
154 // ----- Get static instance ---------------------------------------------
156  if (!fgInstance) fgInstance = new CbmStsPhysics();
157  return fgInstance;
158 }
159 // -------------------------------------------------------------------------
160 
161 
162 // ----- Interpolate a value from a data table -------------------------
163 Double_t CbmStsPhysics::InterpolateDataTable(Double_t eEquiv,
164  map<Double_t, Double_t>& table) {
165 
166  std::map<Double_t, Double_t>::iterator it = table.lower_bound(eEquiv);
167 
168  // Input value smaller than or equal to first table entry:
169  // return first value
170  if (it == table.begin()) return it->second;
171 
172  // Input value larger than last table entry: return last value
173  if (it == table.end()) return (--it)->second;
174 
175  // Else: interpolate from table values
176  Double_t e2 = it->first;
177  Double_t v2 = it->second;
178  it--;
179  Double_t e1 = it->first;
180  Double_t v1 = it->second;
181  return (v1 + (eEquiv - e1) * (v2 - v1) / (e2 - e1));
182 }
183 // -------------------------------------------------------------------------
184 
185 
186 // ----- Landau Width ------------------------------------------------
187 Double_t CbmStsPhysics::LandauWidth(Double_t mostProbableCharge) {
188 
189  // --- Get interpolated value from the data table
190  return InterpolateDataTable(mostProbableCharge, fLandauWidth);
191 }
192 // -------------------------------------------------------------------------
193 
194 
195 // ----- Particle charge for PDG PID ----------------------------------
196 Double_t CbmStsPhysics::ParticleCharge(Int_t pid) {
197 
198  Double_t charge = 0.;
199 
200  // --- For particles in the TDatabasePDG. Note that TParticlePDG
201  // --- gives the charge in units of |e|/3, God knows why.
202  TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pid);
203  if (particle) charge = particle->Charge() / 3.;
204 
205  // --- For ions
206  else if (pid > 1000000000 && pid < 1010000000) {
207  Int_t myPid = pid / 10000;
208  charge = Double_t(myPid - (myPid / 1000) * 1000);
209  }
210 
211  return charge;
212 }
213 // -------------------------------------------------------------------------
214 
215 
216 // ----- Particle mass for PDG PID ------------------------------------
217 Double_t CbmStsPhysics::ParticleMass(Int_t pid) {
218 
219  Double_t mass = -1.;
220 
221  // --- For particles in the TDatabasePDG
222  TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pid);
223  if (particle) mass = particle->Mass();
224 
225  // --- For ions
226  else if (pid > 1000000000 && pid < 1010000000) {
227  Int_t myPid = pid - 1e9;
228  myPid -= (myPid / 10000) * 10000;
229  mass = Double_t(myPid / 10);
230  }
231 
232  return mass;
233 }
234 // -------------------------------------------------------------------------
235 
236 
237 // ----- Read data tables for stopping power ---------------------------
239 
240  // The table with errors for Landau distribution:
241  // MP charge (e) --> half width of charge distribution (e)
242 
243  TString dir = gSystem->Getenv("VMCWORKDIR");
244  TString errFileName = dir + "/parameters/sts/LandauWidthTable.txt";
245 
246  ifstream inFile;
247  Double_t q, err;
248 
249  // --- Read electron stopping power
250  inFile.open(errFileName);
251  if (inFile.is_open()) {
252  while (true) {
253  inFile >> q;
254  inFile >> err;
255  if (inFile.eof()) break;
256  fLandauWidth[q] = err;
257  }
258  inFile.close();
259  LOG(info) << "StsPhysics: " << setw(5) << right << fLandauWidth.size()
260  << " values read from " << errFileName;
261  } else
262  LOG(fatal) << "StsPhysics: Could not read from " << errFileName;
263 }
264 // -------------------------------------------------------------------------
265 
266 
267 // ----- Read data tables for stopping power ---------------------------
269 
270  // The data tables are obtained from the NIST ESTAR and PSTAR databases:
271  // http://www.nist.gov/pml/data/star/index.cfm
272  // The first column in the tables is the kinetic energy in MeV, the second
273  // one is the specific stopping power in MeV*cm^2/g for Silicon.
274  // Internally, the values are stored in GeV and GeV*cm^2/g, respectively.
275  // Conversion MeV->GeV is done when reading from file.
276 
277  TString dir = gSystem->Getenv("VMCWORKDIR");
278  TString eFileName = dir + "/parameters/sts/dEdx_Si_e.txt";
279  TString pFileName = dir + "/parameters/sts/dEdx_Si_p.txt";
280 
281  ifstream inFile;
282  Double_t e, dedx;
283 
284  // --- Read electron stopping power
285  inFile.open(eFileName);
286  if (inFile.is_open()) {
287  while (true) {
288  inFile >> e;
289  inFile >> dedx;
290  if (inFile.eof()) break;
291  e *= 1.e-3; // MeV -> GeV
292  dedx *= 1.e-3; // MeV -> GeV
293  fStoppingElectron[e] = dedx;
294  }
295  inFile.close();
296  LOG(info) << "StsPhysics: " << setw(5) << right << fStoppingElectron.size()
297  << " values read from " << eFileName;
298  } else
299  LOG(fatal) << "StsPhysics: Could not read from " << eFileName;
300 
301  // --- Read proton stopping power
302  inFile.open(pFileName);
303  if (inFile.is_open()) {
304  while (true) {
305  inFile >> e;
306  inFile >> dedx;
307  if (inFile.eof()) break;
308  e *= 1.e-3; // MeV -> GeV
309  dedx *= 1.e-3; // MeV -> GeV
310  fStoppingProton[e] = dedx;
311  }
312  inFile.close();
313  LOG(info) << "StsPhysics: " << setw(5) << right << fStoppingProton.size()
314  << " values read from " << pFileName;
315  } else
316  LOG(fatal) << "StsPhysics: Could not read from " << pFileName;
317 }
318 // -------------------------------------------------------------------------
319 
320 
321 // ----- Set the parameters for the Urban model ------------------------
323 
324  // --- Mean ionisation potential according to PHYS333 2.1
325  fUrbanI = 1.6e-8 * TMath::Power(z, 0.9); // in GeV
326 
327  // --- Maximal energy loss (delta electron threshold) [GeV]
328  // TODO: 1 MeV is the default setting in our transport simulation.
329  // It would be desirable to obtain the actually used value automatically.
330  fUrbanEmax = 1.e-3;
331 
332  // --- The following parameters were defined according the GEANT3 choice
333  // --- described in PHYS332 2.4 (p.264)
334 
335  // --- Oscillator strengths of energy levels
336  fUrbanF1 = 1. - 2. / z;
337  fUrbanF2 = 2. / z;
338 
339  // --- Energy levels [GeV]
340  fUrbanE2 = 1.e-8 * z * z;
341  fUrbanE1 =
342  TMath::Power(fUrbanI / TMath::Power(fUrbanE2, fUrbanF2), 1. / fUrbanF1);
343 
344  // --- Relative weight excitation / ionisation
345  fUrbanR = 0.4;
346 
347  // --- Screen output
348  LOG(info) << "StsPhysics: Urban parameters for z = " << z << " :";
349  LOG(info) << "I = " << fUrbanI * 1.e9 << " eV, Emax = " << fUrbanEmax * 1.e9
350  << " eV, E1 = " << fUrbanE1 * 1.e9
351  << " eV, E2 = " << fUrbanE2 * 1.e9 << " eV, f1 = " << fUrbanF1
352  << ", f2 = " << fUrbanF2 << ", r = " << fUrbanR;
353 }
354 // -------------------------------------------------------------------------
355 
356 
357 // ----- Stopping power ------------------------------------------------
358 Double_t CbmStsPhysics::StoppingPower(Double_t eKin, Int_t pid) {
359 
360  Double_t mass = ParticleMass(pid);
361  if (mass < 0.) return 0.;
362  Double_t charge = ParticleCharge(pid);
363  Bool_t isElectron = (pid == 11 || pid == -11);
364 
365  return StoppingPower(eKin, mass, charge, isElectron);
366 }
367 // -------------------------------------------------------------------------
368 
369 
370 // ----- Stopping power ------------------------------------------------
371 Double_t CbmStsPhysics::StoppingPower(Double_t energy,
372  Double_t mass,
373  Double_t charge,
374  Bool_t isElectron) {
375 
376  // --- Get interpolated value from data table
377  Double_t stopPower = -1.;
378  if (isElectron)
379  stopPower = InterpolateDataTable(energy, fStoppingElectron);
380  else {
381  Double_t eEquiv = energy * kProtonMass / mass; // equiv. proton energy
382  stopPower = InterpolateDataTable(eEquiv, fStoppingProton);
383  }
384 
385  // --- Calculate stopping power (from specific SP and density of silicon)
386  stopPower *= kSiDensity; // density of silicon
387 
388  // --- For non-electrons: scale with squared charge
389  if (!isElectron) stopPower *= (charge * charge);
390 
391  return stopPower;
392 }
393 // -------------------------------------------------------------------------
CbmStsPhysics
Auxiliary class for physics processes in Silicon.
Definition: CbmStsPhysics.h:24
CbmStsPhysics::ReadDataTablesStoppingPower
void ReadDataTablesStoppingPower()
Read stopping power data table from file.
Definition: CbmStsPhysics.cxx:268
CbmStsPhysics::ParticleCharge
static Double_t ParticleCharge(Int_t pid)
Particle charge from PDG particle ID.
Definition: CbmStsPhysics.cxx:196
CbmSts
Definition: CbmStsDefs.h:15
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmStsPhysics::LandauWidth
Double_t LandauWidth(Double_t mostProbableCharge)
Half width at half max of Landau distribution in ultra-relativistic case.
Definition: CbmStsPhysics.cxx:187
CbmStsPhysics::~CbmStsPhysics
virtual ~CbmStsPhysics()
Definition: CbmStsPhysics.cxx:53
CbmStsPhysics::fgInstance
static CbmStsPhysics * fgInstance
Singleton instance.
Definition: CbmStsPhysics.h:194
CbmSts::kSiCharge
const Double_t kSiCharge
Silicon charge [e].
Definition: CbmStsDefs.h:18
CbmStsPhysics::CbmStsPhysics
CbmStsPhysics()
Constructor.
CbmStsPhysics::StoppingPower
Double_t StoppingPower(Double_t eKin, Int_t pid)
Stopping power (average specific energy loss) in Silicon.
Definition: CbmStsPhysics.cxx:358
CbmStsPhysics::Instance
static CbmStsPhysics * Instance()
Accessor to singleton instance.
Definition: CbmStsPhysics.cxx:155
CbmStsDefs.h
CbmStsPhysics::ParticleMass
static Double_t ParticleMass(Int_t pid)
Particle mass from PDG particle ID.
Definition: CbmStsPhysics.cxx:217
d
double d
Definition: P4_F64vec2.h:24
log
friend F32vec4 log(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:135
CbmStsPhysics::EnergyLoss
Double_t EnergyLoss(Double_t dz, Double_t mass, Double_t eKin, Double_t dedx) const
Energy loss in a Silicon layer.
Definition: CbmStsPhysics.cxx:112
CbmStsPhysics.h
ClassImp
ClassImp(CbmConverterManager) InitStatus CbmConverterManager
Definition: CbmConverterManager.cxx:12
CbmStsPhysics::ReadDataTablesLandauWidth
void ReadDataTablesLandauWidth()
Read Landau width data table from file.
Definition: CbmStsPhysics.cxx:238
CbmStsPhysics::SetUrbanParameters
void SetUrbanParameters(Double_t z)
Calculate the parameters for the Urban model.
Definition: CbmStsPhysics.cxx:322
CbmSts::kSiDensity
const Double_t kSiDensity
Silicon density [g/cm^3].
Definition: CbmStsDefs.h:22
CbmStsPhysics::InterpolateDataTable
Double_t InterpolateDataTable(Double_t eKin, std::map< Double_t, Double_t > &table)
Interpolate a value from the data tables.
Definition: CbmStsPhysics.cxx:163
CbmStsPhysics::ElectricField
static Double_t ElectricField(Double_t vBias, Double_t vFd, Double_t dZ, Double_t z)
Electric field magnitude in a silicon sensor as function of z.
Definition: CbmStsPhysics.cxx:102
CbmStsPhysics::DiffusionWidth
static Double_t DiffusionWidth(Double_t z, Double_t d, Double_t vBias, Double_t vFd, Double_t temperature, Int_t chargeType)
Definition: CbmStsPhysics.cxx:58
CbmSts::kProtonMass
const Double_t kProtonMass
Proton mass [GeV].
Definition: CbmStsDefs.h:26