feat: current simulation

This commit is contained in:
liuyihui 2023-06-01 15:16:26 +08:00
parent 9fac823a78
commit 71e6945cef
5 changed files with 333 additions and 138 deletions

View File

@ -23,23 +23,24 @@ int main(int argc, char* argv[]) {
TApplication app("app", &argc, argv);
MediumMagboltz gas;
// gas.LoadGasFile("helium3_100atm_3000K.gas");
// gas.LoadIonMobility("IonMobility_He+_He.txt");
gas.LoadGasFile("helium3_100atm_3000K.gas");
gas.LoadIonMobility("IonMobility_He+_He.txt");
gas.SetComposition("helium-3");
gas.SetPressure(10.0 * 760.);
gas.SetTemperature(300);
gas.SetFieldGrid(100, 15e4, 20, true);
gas.SetFieldGrid(1e5, 1.35e5, 3, true);
gas.EnableThermalMotion();
gas.GenerateGasTable(10);
gas.WriteGasFile("helium3_100atm_3000K.gas");
// gas.WriteGasFile("helium3_100atm_3000K.gas");
TCanvas* cM = new TCanvas("cM", "", 600, 600);
ViewMedium mediumView;
mediumView.SetMedium(&gas);
mediumView.SetCanvas(cM);
mediumView.PlotElectronVelocity('e');
mediumView.PlotHoleVelocity('e', true);
// mediumView.PlotIonVelocity();
// mediumView.PlotHoleVelocity('e', true);
app.Run(kTRUE);
return 0;
@ -54,15 +55,15 @@ int main(int argc, char* argv[]) {
cmp.AddTube(rTube, vTube, 0, "t");
cmp.AddReadout("s");
// ViewField fieldView;
// fieldView.SetComponent(&cmp);
// fieldView.SetPlane(0, 0, 1, 0, 0, 0);
// fieldView.SetArea(-0.5, -0.5, 0.5, 0.5);
// fieldView.SetVoltageRange(0., 1000.);
// fieldView.PlotContour("e");
ViewField fieldView;
fieldView.SetComponent(&cmp);
fieldView.SetPlane(0, 0, 1, 0, 0, 0);
fieldView.SetArea(-0.5, -0.5, 0.5, 0.5);
fieldView.SetVoltageRange(0., 1000.);
fieldView.PlotContour("e");
// fieldView.PlotProfile(rWire, 0, 0, rTube, 0, 0, "e");
// app.Run(kTRUE);
// return 0;
app.Run(kTRUE);
return 0;
const double tStep = 0.5;
Sensor sensor;

57
mdt.cpp
View File

@ -16,7 +16,7 @@
#include <fstream>
#include <iostream>
#define f_PLOT 0
#define f_PLOT 1
#define f_PRINT 2
using namespace Garfield;
@ -28,7 +28,11 @@ const double vTube = 0.;
const double tStep = 1;
const int nbins = 10000;
const int N = 1000;
const int N = 1;
#if f_PRINT
std::ofstream outfile;
#endif
bool readTransferFunction(Sensor& sensor, std::string filename) {
std::ifstream infile;
@ -73,11 +77,12 @@ double riseTime(Sensor& sensor) {
return 1.e-9 * (tr - tl) * tStep;
}
void randPos(double* x, double* y) {
void randPos(double* x, double* y, double* z) {
double a1 = 2 * M_PI * (rand() / (double)RAND_MAX);
double a2 = (rand() / (double)RAND_MAX) * rTube * 0.95;
double a2 = (rand() / (double)RAND_MAX) * rTube;
*x = cos(a1) * a2 + rWire;
*y = sin(a1) * a2 + rWire;
*z = (2 * (rand() / (double)RAND_MAX) - 1) * 10;
}
void randDir(double* x, double* y, double* z) {
@ -101,12 +106,19 @@ std::tuple<int, double> driftElectron(TrackSrim track, DriftLineRKF drift) {
drift.DriftIon(xc, yc, zc, tc);
drift.SetElectronSignalScalingFactor(nc);
drift.DriftElectron(xc, yc, zc, tc);
#if f_PRINT
outfile << xc << " " << yc << " " << zc << " " << tc << std::endl;
#endif
}
return {ne, eDep};
}
int main(int argc, char* argv[]) {
#if f_PRINT
outfile.open("test.txt", std::ios::out);
#endif
srand(time(NULL));
TApplication app("app", &argc, argv);
@ -162,52 +174,51 @@ int main(int argc, char* argv[]) {
alpha.EnablePlotting(&driftView);
#endif
#if f_PRINT
std::ofstream outfile;
outfile.open("proton.txt", std::ios::out);
#endif
double ek, rt, px, py, dx, dy, dz;
double ek, rt, px, py, pz, dx, dy, dz;
for (int i = 0; i < N; i++) {
#if f_PLOT
driftView.Clear();
#endif
sensor.ClearSignal();
randPos(&px, &py);
randPos(&px, &py, &pz);
randDir(&dx, &dy, &dz);
ek = 0.6e6 * (rand() / (double)RAND_MAX) + 0.2e6;
ek = 0.8e6 * (rand() / (double)RAND_MAX) + 0.2e6;
// proton.SetKineticEnergy(ek);
proton.NewTrack(px, py, pz, 0, dx, dy, dz);
// tritium.SetKineticEnergy(ek);
tritium.NewTrack(px, py, pz, 0, -dx, -dy, -dz);
// alpha.SetKineticEnergy(ek);
// alpha.NewTrack(px, py, pz, 0, dx, dy, dz);
proton.SetKineticEnergy(ek);
proton.NewTrack(px, py, 0, 0, dx, dy, dz);
auto [ne, eDep] = driftElectron(proton, drift);
// tritium.NewTrack(px, py, 0, 0, -dx, -dy, -dz);
// alpha.NewTrack(rTrack, rTrack, 0, 0, -1, 1, 0);
auto [ne2, eDep2] = driftElectron(tritium, drift);
// sensor.IntegrateSignals();
sensor.ConvoluteSignals();
rt = riseTime(sensor) * 1e9;
#if f_PLOT
// cD->Clear();
cD->Clear();
cellView.Plot2d();
driftView.Plot(true, false);
cD->Update();
signalView.SetSensor(&sensor);
// signalView.PlotSignal("s", "tei");
// signalView.PlotSignal("s", "ei");
signalView.PlotSignal("s", "t");
#endif
#if f_PRINT > 1
printf("Ek = %.3f, Rise Time = %.3f\n", ek, rt);
printf("Energy Deposit = %.3f, Rise Time = %.3f\n", eDep + eDep2, rt);
printf("Position = (%.2f, %.2f), Direction = (%.2f, %.2f, %.2f)\n", px, py, dx, dy, dz);
printf("Total Count = %d, Energy Deposit = %.3f, W = %.3f\n", ne, eDep, eDep / ne);
#endif
#if f_PRINT
outfile << eDep << " " << rt << std::endl;
// outfile << eDep + eDep2 << " " << rt << std::endl;
outfile << sensor.GetSignal("s", 0);
for (int i = 1; i < nbins; i++) outfile << " " << sensor.GetSignal("s", i);
outfile << std::endl;
#endif
}

View File

@ -1,100 +1,100 @@
0.000000 0.000000
0.008333 0.039966
0.016667 0.076670
0.025000 0.110312
0.033333 0.141080
0.041667 0.169153
0.050000 0.194700
0.058333 0.217880
0.066667 0.238844
0.075000 0.257733
0.083333 0.274684
0.091667 0.289821
0.100000 0.303265
0.108333 0.315130
0.116667 0.325521
0.125000 0.334538
0.133333 0.342278
0.141667 0.348829
0.150000 0.354275
0.158333 0.358695
0.166667 0.362165
0.175000 0.364754
0.183333 0.366529
0.191667 0.367551
0.200000 0.367879
0.200000 0.367879
0.222222 0.365770
0.244444 0.360036
0.266667 0.351463
0.288889 0.340711
0.311111 0.328334
0.333333 0.314793
0.355556 0.300468
0.377778 0.285675
0.400000 0.270671
0.400000 0.270671
0.416667 0.259405
0.433333 0.248211
0.450000 0.237148
0.466667 0.226268
0.483333 0.215611
0.500000 0.205212
0.516667 0.195098
0.533333 0.185289
0.550000 0.175802
0.566667 0.166647
0.583333 0.157832
0.600000 0.149361
0.616667 0.141236
0.633333 0.133456
0.650000 0.126016
0.666667 0.118913
0.683333 0.112141
0.700000 0.105691
0.716667 0.099556
0.733333 0.093726
0.750000 0.088192
0.766667 0.082943
0.783333 0.077970
0.800000 0.073263
0.800000 0.073263
0.830769 0.065232
0.861538 0.058001
0.892308 0.051507
0.923077 0.045685
0.953846 0.040476
0.984615 0.035824
1.015385 0.031675
1.046154 0.027982
1.076923 0.024697
1.107692 0.021780
1.138462 0.019193
1.169231 0.016901
1.200000 0.014873
1.230769 0.013079
1.261538 0.011494
1.292308 0.010095
1.323077 0.008862
1.353846 0.007775
1.384615 0.006818
1.415385 0.005976
1.446154 0.005235
1.476923 0.004584
1.507692 0.004012
1.538462 0.003510
1.569231 0.003070
1.600000 0.002684
1.630769 0.002345
1.661538 0.002049
1.692308 0.001789
1.723077 0.001562
1.753846 0.001363
1.784615 0.001189
1.815385 0.001037
1.846154 0.000904
1.876923 0.000788
1.907692 0.000687
1.938462 0.000599
1.969231 0.000521
2.000000 0.000454
0.004167 0.039966
0.008333 0.076670
0.012500 0.110312
0.016667 0.141080
0.020833 0.169153
0.025000 0.194700
0.029167 0.217880
0.033333 0.238844
0.037500 0.257733
0.041667 0.274684
0.045833 0.289821
0.050000 0.303265
0.054167 0.315130
0.058333 0.325521
0.062500 0.334538
0.066667 0.342278
0.070833 0.348829
0.075000 0.354275
0.079167 0.358695
0.083333 0.362165
0.087500 0.364754
0.091667 0.366529
0.095833 0.367551
0.100000 0.367879
0.100000 0.367879
0.111111 0.365770
0.122222 0.360036
0.133333 0.351463
0.144444 0.340711
0.155556 0.328334
0.166667 0.314793
0.177778 0.300468
0.188889 0.285675
0.200000 0.270671
0.200000 0.270671
0.208333 0.259405
0.216667 0.248211
0.225000 0.237148
0.233333 0.226268
0.241667 0.215611
0.250000 0.205212
0.258333 0.195098
0.266667 0.185289
0.275000 0.175802
0.283333 0.166647
0.291667 0.157832
0.300000 0.149361
0.308333 0.141236
0.316667 0.133456
0.325000 0.126016
0.333333 0.118913
0.341667 0.112141
0.350000 0.105691
0.358333 0.099556
0.366667 0.093726
0.375000 0.088192
0.383333 0.082943
0.391667 0.077970
0.400000 0.073263
0.400000 0.073263
0.415385 0.065232
0.430769 0.058001
0.446154 0.051507
0.461538 0.045685
0.476923 0.040476
0.492308 0.035824
0.507692 0.031675
0.523077 0.027982
0.538462 0.024697
0.553846 0.021780
0.569231 0.019193
0.584615 0.016901
0.600000 0.014873
0.615385 0.013079
0.630769 0.011494
0.646154 0.010095
0.661538 0.008862
0.676923 0.007775
0.692308 0.006818
0.707692 0.005976
0.723077 0.005235
0.738462 0.004584
0.753846 0.004012
0.769231 0.003510
0.784615 0.003070
0.800000 0.002684
0.815385 0.002345
0.830769 0.002049
0.846154 0.001789
0.861538 0.001562
0.876923 0.001363
0.892308 0.001189
0.907692 0.001037
0.923077 0.000904
0.938462 0.000788
0.953846 0.000687
0.969231 0.000599
0.984615 0.000521
1.000000 0.000454

View File

@ -7,15 +7,15 @@ def step(x):
Q = 1e-15 # 1 fC
Cf = 1e-12 # 1 pF
tau = 2e-7 # 0.2 us
tau = 1e-7 # 0.1 us
t = np.append(np.linspace(0, tau, 25), np.linspace(tau, 2 * tau, 10))
t = np.append(t, np.linspace(2 * tau, 4 * tau, 25))
t = np.append(t, np.linspace(4 * tau, 10 * tau, 40))
v = Q / Cf * (t / tau) * np.exp(-t / tau) * step(t)
# plt.plot(t, v, 'o', markersize=2)
# plt.xlabel("time")
# plt.plot(t * 1e6, v, 'o', markersize=2)
# plt.xlabel("time [us]")
# plt.ylabel("amplitude")
# plt.show()

183
test.txt Normal file

File diff suppressed because one or more lines are too long