10 #include <FairLogger.h>
12 #include <TDatabasePDG.h>
14 #include <TParticlePDG.h>
42 LOG(info) <<
"Instantiating STS Physics... ";
43 ReadDataTablesStoppingPower();
44 ReadDataTablesLandauWidth();
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 <<
")";
74 if (temperature < 0.) {
75 LOG(error) <<
"StsPhysics: illegal temperature value " << temperature;
81 Double_t diffConst = 8.61733e-5 * temperature;
86 if (chargeType == 0) {
87 tau = 0.5 *
d *
d / vFd
88 *
log((vBias + (1. - 2. * z /
d) * vFd) / (vBias - vFd));
89 }
else if (chargeType == 1) {
90 tau = -0.5 *
d *
d / vFd *
log(1. - 2. * vFd * z /
d / (vBias + vFd));
92 LOG(error) <<
"StsPhysics: Illegal charge type " << chargeType;
96 return sqrt(2. * diffConst * tau);
106 return (vBias + vFd * (2. * z / dZ - 1.)) / dZ;
115 Double_t dedx)
const {
118 Double_t gamma = (eKin + mass) / mass;
119 Double_t beta2 = 1. - 1. / (gamma * gamma);
122 Double_t xAux = 2. * mass * beta2 * gamma * gamma;
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);
137 Int_t n1 = gRandom->Poisson(sigma1 * dz);
138 Int_t n2 = gRandom->Poisson(sigma2 * dz);
139 Int_t n3 = gRandom->Poisson(sigma3 * dz);
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));
149 return (n1 * fUrbanE1 + n2 * fUrbanE2 + eLossIon);
164 map<Double_t, Double_t>& table) {
166 std::map<Double_t, Double_t>::iterator it = table.lower_bound(eEquiv);
170 if (it == table.begin())
return it->second;
173 if (it == table.end())
return (--it)->second;
176 Double_t e2 = it->first;
177 Double_t v2 = it->second;
179 Double_t e1 = it->first;
180 Double_t v1 = it->second;
181 return (v1 + (eEquiv - e1) * (v2 - v1) / (e2 - e1));
190 return InterpolateDataTable(mostProbableCharge, fLandauWidth);
198 Double_t charge = 0.;
202 TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pid);
203 if (particle) charge = particle->Charge() / 3.;
206 else if (pid > 1000000000 && pid < 1010000000) {
207 Int_t myPid = pid / 10000;
208 charge = Double_t(myPid - (myPid / 1000) * 1000);
222 TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pid);
223 if (particle) mass = particle->Mass();
226 else if (pid > 1000000000 && pid < 1010000000) {
227 Int_t myPid = pid - 1e9;
228 myPid -= (myPid / 10000) * 10000;
229 mass = Double_t(myPid / 10);
243 TString dir = gSystem->Getenv(
"VMCWORKDIR");
244 TString errFileName = dir +
"/parameters/sts/LandauWidthTable.txt";
250 inFile.open(errFileName);
251 if (inFile.is_open()) {
255 if (inFile.eof())
break;
256 fLandauWidth[q] = err;
259 LOG(info) <<
"StsPhysics: " << setw(5) << right << fLandauWidth.size()
260 <<
" values read from " << errFileName;
262 LOG(fatal) <<
"StsPhysics: Could not read from " << errFileName;
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";
285 inFile.open(eFileName);
286 if (inFile.is_open()) {
290 if (inFile.eof())
break;
293 fStoppingElectron[e] = dedx;
296 LOG(info) <<
"StsPhysics: " << setw(5) << right << fStoppingElectron.size()
297 <<
" values read from " << eFileName;
299 LOG(fatal) <<
"StsPhysics: Could not read from " << eFileName;
302 inFile.open(pFileName);
303 if (inFile.is_open()) {
307 if (inFile.eof())
break;
310 fStoppingProton[e] = dedx;
313 LOG(info) <<
"StsPhysics: " << setw(5) << right << fStoppingProton.size()
314 <<
" values read from " << pFileName;
316 LOG(fatal) <<
"StsPhysics: Could not read from " << pFileName;
325 fUrbanI = 1.6e-8 * TMath::Power(z, 0.9);
336 fUrbanF1 = 1. - 2. / z;
340 fUrbanE2 = 1.e-8 * z * z;
342 TMath::Power(fUrbanI / TMath::Power(fUrbanE2, fUrbanF2), 1. / fUrbanF1);
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;
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);
365 return StoppingPower(eKin, mass, charge, isElectron);
377 Double_t stopPower = -1.;
379 stopPower = InterpolateDataTable(energy, fStoppingElectron);
382 stopPower = InterpolateDataTable(eEquiv, fStoppingProton);
389 if (!isElectron) stopPower *= (charge * charge);