Hi,
I have reproduce my programme, but the result is still the same as below
would you please help me to find out the reason, Thanks a lot.
would you please help me to find out the reason, Thanks a lot.
You need to remove the lines
drift.SetHoleSignalScalingFactor(ne);
drift.DriftHole(xc, yc, zc, tc);
or (better) replace Hole by Ion. The TRIM example in Examples/Srim uses silicon but your program uses gas as detection medium.
To simulate ion drift lines, you will also need to import a table of ion mobilities (see e. g. the drift tube example or the GEM example).
Hi,
Thanks a lot for your reply
I don’t need the ion drift line, I only want to get the induced current, I will try to delete.
Hi,
This is my code
//induced current
#include <iostream>
#include <TApplication.h>
#include <TH1F.h>
#include <TCanvas.h>
#include "Garfield/MediumSilicon.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackTrim.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/MediumMagboltz.hh"
using namespace Garfield;
int main(int argc, char *argv[]) {
// Application
TApplication app("app", &argc, argv);
ComponentAnsys123 fm;
fm.Initialise("ELIST.lis", "NLIST.lis", "MPLIST.lis", "PRNSOL.lis", "cm");
fm.SetWeightingField("weight1.lis", "readout");
double x0, y0, z0, x1, y1, z1;
fm.GetBoundingBox(x0, y0, z0, x1, y1, z1);
// Define the medium.
MediumMagboltz gas;
gas.LoadGasFile("Ar_90_CH4_10.gas");
fm.SetMedium(0, &gas);
Sensor sensor;
sensor.AddComponent(&fm);
sensor.AddElectrode(&fm, "readout");
// Set the time bins for the induced current.
const unsigned int nTimeBins = 9000;
const double tmin = 1000.;
const double tmax = 100000.;
const double tstep = (tmax - tmin) / nTimeBins;
sensor.SetTimeWindow(tmin, tstep, nTimeBins);
// Read the TRIM output file.
TrackTrim tr;
const std::string filename = "EXYZ-P10.txt";
// Import the first 100 ions.
if (!tr.ReadFile(filename, 100)) {
std::cerr << "Reading TRIM EXYZ file failed.\n";
return 1;
}
tr.SetWorkFunction(26);
tr.SetFanoFactor(0.179);
// Connect the track to a sensor.
tr.SetSensor(&sensor);
ViewSignal signalView;
signalView.SetSensor(&sensor);
signalView.PlotSignal("readout", true, true, true);
app.Run();
return 0;
}
Apologies for the delay in getting back to you.
I just tried to load the weight1.lis you sent by e-mail but it’s not compatible with the latest set of ANSYS field map files (ELIST.lis, NLIST.lis, PRNSOL.lis, etc.) I have from you.
Can you try to do a bit of debugging to identify the cause of the problem:
- Make a plot of the weighting potential. Does it look like you expect it to?
- Simulate a few electron and/or ion drift lines along a track. Do either the electrons or the ions go to the readout electrode?
Also as I tried to suggest in this thread, if the signal you measure is proportional to the collected charge, you could simplify the calculation by simply counting the number of electrons/ions arriving at each electrode.
1 Like
HI,
Thanks for your reply
That’s really kind of you for helping me to solve lots of problems.And without that I may finally give up.
Yep,I finally found that the four lis files is not compatible with the weight.lis file, So I reproduce all these files , but the result is still the same like this in the below.
Besides, the thickness of my detector is 30mm at least, but after importing it to the Garfield++, it turns to be 0.3mm, however the potential plot seems right.That’s wired. Is that the reason for the wrong units setting? I have tried to simplified the model to three parts ,one part is a metal , and along the up and down of the metal , there are two gas areas, the metal is used to be the electrode to collect the electron. And then I will set the metal’s two areas to 1V, and the other all area about gas to 0 V, and this will be the new weight.lis. I don’t know whether it can work.
Besides, how to get the number of electrons/ions arriving at each electrode. Which code or output files can make it
- There might be an issue with the weighting field map or a bug in the code but to look into it I would need a consistent set of Ansys output files (
.lisfiles). - The Ansys output files do not include the unit in which the mesh coordinates are given. So you need to specify the length unit when you import the files in Garfield++ using the
Initialisefunction (it’s the last argument of the function, in the code you posted earlier it was"cm").
Sorry, I realise just now that in this program you don’t simulate any drift lines.
You need to simulate electron or ion drift lines to calculate the induced signal.
In your program you have a line like this to simulate an electron drift line:
drift.DriftElectron(xc, yc, zc, tc);
You can get the end point of this drift line (i. e. the point at which the electron stopped drifting) using something like this
double xf, yf, zf, tf;
int status = 0;
drift.GetEndPoint(xf, yf, zf, tf, status);
Now, knowing the coordinates of the end point xf, yf, zf, you can determine if it was collected by the electrode that you want to read out or not. For instance, if the readout electrode is the plane y = 0 (but this is really just an example, you’ll need to adapt it):
if (fabs(yf) <= 1.e-6) {
// Electron ended up on the readout plane (within some tolerance).
// Increase a counter.
++nElectronsOnReadoutPlane;
// ...
}
Another useful function is GetAvalancheSize which will tell you how many electrons and ions were produced in the avalanche (if there is multiplication in your detector). Extending the example above:
drift.DriftElectron(xc, yc, zc, tc);
double xf, yf, zf, tf;
int status = 0;
drift.GetEndPoint(xf, yf, zf, tf, status);
double neaval = 0., niaval = 0.;
drift.GetAvalancheSize(neaval, niaval);
if (fabs(yf) <= 1.e-6) {
// Electron ended up on the readout plane (within some tolerance).
// Increase a counter.
nElectronsOnReadoutPlane += neval;
// ...
}
Ok, I konw it now. I have simulated the electron and drift lines before, but there is only one or two clusters generated, and the drift lines only has one or two. I think maybe it is the reason for the low pressure. But now I think the reason is from the gas file generated, the min and max value of electric field is for one Torr only. I will try to solve it and reproduce the programme.
Ok, so maybe I should set the length unit to mm
Ah, I realize it now. If I can get the induced current, then I would not need the number of them. The code below is what I used to get the drift lines. With these drift lines I can then get the induced current ,am I right?
*#include <iostream>
#include <cstdlib>
#include <TH1F.h>
#include <TROOT.h>
#include <TApplication.h>
#include <TCanvas.h>
#include "Garfield/FundamentalConstants.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewFEMesh.hh"
#include "Garfield/ViewMedium.hh"
#include "Garfield/TrackSrim.hh"
#include "Garfield/Random.hh"
#include "Garfield/Plotting.hh"
#include "Garfield/ComponentConstant.hh"
using namespace Garfield;
int main(int argc, char * argv[]) {
TApplication app("app", &argc, argv);
ComponentAnsys123 fm;
fm.Initialise("ELIST.lis", "NLIST.lis", "MPLIST.lis", "PRNSOL.lis", "cm");
fm.PrintRange();
fm.PrintMaterials();
// Get the extent of the field map.
double x0, y0, z0, x1, y1, z1;
fm.GetBoundingBox(x0, y0, z0, x1, y1, z1);
ViewField fieldView;
constexpr bool plotField = true;
if (plotField) {
fieldView.SetComponent(&fm);
// Set the normal vector of the viewing plane.
fieldView.SetPlane(0, 0, -1, 0, 0, 0.5 * (z0 + z1));
fieldView.PlotContour();
}
Sensor sensor;
sensor.AddComponent(&fm);
MediumMagboltz gas;
gas.LoadGasFile("4He.gas");
//gas.LoadGasFile("ar_93_co2_7.gas");
fm.SetMedium(0, &gas);
fm.PrintMaterials();
TrackHeed track;
track.SetParticle("proton");
track.SetEnergy(100.e6);
track.SetSensor(&sensor);
track.Initialise(&gas, true);
DriftLineRKF drift;
drift.SetSensor(&sensor);
ViewDrift driftview;
TCanvas* CD = nullptr;
constexpr bool plotDrift = true;
if (plotDrift) {
CD = new TCanvas("CD", "", 600, 600);
driftview.SetCanvas(CD);
drift.EnablePlotting(&driftview);
track.EnablePlotting(&driftview);
tr.EnablePlotting(&driftView);
}
const unsigned int nTracks = 1;
size_t nClusters = 0;
for (unsigned int j = 0; j < nTracks; ++j) {
sensor.ClearSignal();
// Starting point of the track.
double xt = -0.1, yt = 0.204, zt = -0.07, tt = 0.;
// Initial direction of the track.
double dxt = 0., dyt = -1., dzt = 0.;
track.NewTrack(xt, yt, zt, tt, dxt, dyt, dzt);
// Loop over the clusters on the track.
double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.;
int nc = 0;
while (track.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
++nClusters;
// Loop over the electrons in the cluster.
for (int k = 0; k < nc; k++) {
double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.;
double dxe = 0., dye = 0., dze = 0.;
track.GetElectron(k, xe, ye, ze, te, ee, dxe, dye, dze);
drift.DriftElectron(xe, ye, ze, te);
}
}
}
std::cout << nClusters << " clusters...\n";
const bool twod = true;
const bool axis = true;
driftview.Plot(twod, axis);
app.Run(true);
}
No, there are a few things missing:
- You need to load a map of the weighting potential
fm.SetWeightingField("weight1.lis", "readout");
- You need to instruct the “Sensor” to use this weighting potential/field for computing the induced current:
sensor.AddElectrode(&fm, "readout");
- Typically, you also want to set the time window (
SetTimeWindow) for the signal calculation to make sure it’s in a reasonable range and has a reasonable resolution.
You already had these steps in a previous version of your program, I don’t know why you deleted them.
Also I noticed that you are using TrackHeed for simulating the ionisation by a charged particle track. That’s ok, but I thought you were using TRIM?
Hi,
Thank you for your reply
Sorry for misunderstanding you, this code is an old version for simulating the drift line, I quote it just to show that I have simulated them before. And the whole version is like this, I don’t know whether it can work, and I am trying now.
#include <iostream>
#include <TApplication.h>
#include <TH1F.h>
#include <TCanvas.h>
#include "Garfield/MediumSilicon.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackTrim.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/MediumMagboltz.hh"
using namespace Garfield;
int main(int argc, char *argv[]) {
// Application
TApplication app("app", &argc, argv);
ComponentAnsys123 fm;
fm.Initialise("ELIST.lis", "NLIST.lis", "MPLIST.lis", "PRNSOL.lis", "cm");
fm.SetWeightingField("weight.lis", "readout");
double x0, y0, z0, x1, y1, z1;
fm.GetBoundingBox(x0, y0, z0, x1, y1, z1);
// Define the medium.
MediumMagboltz gas;
gas.LoadGasFile("4He.gas");
fm.SetMedium(0, &gas);
Sensor sensor;
sensor.AddComponent(&fm);
sensor.AddElectrode(&fm, "readout");
//sensor.SetTimeWindow(tmin, tstep, nTimeBins);
sensor.SetTimeWindow(0., 1000., 100);
sensor.GetSignal("readout",1000);
DriftLineRKF drift;
drift.SetSensor(&sensor);
ViewDrift driftview;
TCanvas* CD = nullptr;
constexpr bool plotDrift = true;
if (plotDrift) {
CD = new TCanvas("CD", "", 600, 600);
driftview.SetCanvas(CD);
drift.EnablePlotting(&driftview);
track.EnablePlotting(&driftview);
tr.EnablePlotting(&driftView);
}
const unsigned int nTracks = 1;
size_t nClusters = 0;
for (unsigned int j = 0; j < nTracks; ++j) {
sensor.ClearSignal();
// Starting point of the track.
double xt = -0.1, yt = 0.204, zt = -0.07, tt = 0.;
// Initial direction of the track.
double dxt = 0., dyt = -1., dzt = 0.;
track.NewTrack(xt, yt, zt, tt, dxt, dyt, dzt);
// Loop over the clusters on the track.
double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.;
int nc = 0;
while (track.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
++nClusters;
// Loop over the electrons in the cluster.
for (int k = 0; k < nc; k++) {
double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.;
double dxe = 0., dye = 0., dze = 0.;
track.GetElectron(k, xe, ye, ze, te, ee, dxe, dye, dze);
drift.DriftElectron(xe, ye, ze, te);
}
}
}
std::cout << nClusters << " clusters...\n";
const bool twod = true;
const bool axis = true;
driftview.Plot(twod, axis);
// Read the TRIM output file.
TrackTrim tr;
const std::string filename = "EXYZ-He.txt";
// Import the first 100 ions.
if (!tr.ReadFile(filename, 500)) {
std::cerr << "Reading TRIM EXYZ file failed.\n";
return 1;
}
tr.SetWorkFunction(26);
tr.SetFanoFactor(0.179);
// Connect the track to a sensor.
tr.SetSensor(&sensor);
ViewSignal signalView;
signalView.SetSensor(&sensor);
signalView.PlotSignal("readout", true, true, true);
app.Run();
return 0;
}
Hi,
I have simplified the model and used the new lis files and weight file to run, the result shows killed. I used the old model, it is still the same noticing killed.
Besides, I found that I have to use the TrackHeed class ,otherwise the code below can’t be used ,and that’s the code for calculating the drift line
const unsigned int nTracks = 1;
size_t nClusters = 0;
for (unsigned int j = 0; j < nTracks; ++j) {
sensor.ClearSignal();
// Starting point of the track.
double xt = -0.1, yt = 0.204, zt = -0.07, tt = 0.;
// Initial direction of the track.
double dxt = 0., dyt = -1., dzt = 0.;
**track.NewTrack(xt, yt, zt, tt, dxt, dyt, dzt);**
// Loop over the clusters on the track.
double xc = 0., yc = 0., zc = 0., tc = 0., ec = 0., extra = 0.;
int nc = 0;
while (**track.GetCluster(xc, yc, zc, tc, nc, ec, extra**)) {
++nClusters;
// Loop over the electrons in the cluster.
for (int k = 0; k < nc; k++) {
double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.;
double dxe = 0., dye = 0., dze = 0.;
**track.GetElectron(k, xe, ye, ze, te, ee, dxe, dye, dze);**
drift.DriftElectron(xe, ye, ze, te);
}
}
}
If I delete the TrackHeed class
the code
while (track.GetCluster(xc, yc, zc, tc, nc, ec, extra)) {
++nClusters;
and the code
track.NewTrack(xt, yt, zt, tt, dxt, dyt, dzt);
will notice error
Yes of course. If you delete the lines of code that create and configure a TrackHeed object, you cannot use this object in your program. I don’t understand what’s the question here.
You can but you don’t have to use TrackHeed if you want to simulate drift lines from a charged-particle track. If you want to use TRIM, you can do something like in this example to which I pointed you earlier and which you used in a previous version of your program.
Hi,
Thanks a lot for your reply
I know it now, the drift line come from my TRIM files right? Now I have better understanding about the Garfield++ and their relationship.
#include <iostream>
#include <TApplication.h>
#include <TH1F.h>
#include <TCanvas.h>
#include "Garfield/MediumSilicon.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackTrim.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/TrackHeed.hh"
using namespace Garfield;
int main(int argc, char *argv[]) {
// Application
TApplication app("app", &argc, argv);
ComponentAnsys123 fm;
fm.Initialise("ELIST.lis", "NLIST.lis", "MPLIST.lis", "PRNSOL.lis", "mm");
fm.SetWeightingField("weight.lis", "readout");
double x0, y0, z0, x1, y1, z1;
fm.GetBoundingBox(x0, y0, z0, x1, y1, z1);
ViewField fieldView;
constexpr bool plotField = true;
if (plotField) {
fieldView.SetComponent(&fm);
// Set the normal vector of the viewing plane.
fieldView.SetPlane(0, 0, -1, 0, 0, 0.5 * (z0 + z1));
fieldView.PlotContour();
//fieldView.PlotFieldLines(xf, yf, zf, true, true, kOrange-3);
}
// Define the medium.
MediumMagboltz gas;
gas.LoadGasFile("He-740Torr.gas");
fm.SetMedium(0, &gas);
Sensor sensor;
sensor.AddComponent(&fm);
sensor.AddElectrode(&fm, "readout");
sensor.SetTimeWindow(0., 1000., 100);
sensor.GetSignal("readout",1000);
// Read the TRIM output file.
TrackTrim tr;
const std::string filename = "EXYZ-He.txt";
// Import the first 100 ions.
if (!tr.ReadFile(filename, 500)) {
std::cerr << "Reading TRIM EXYZ file failed.\n";
return 1;
}
tr.SetWorkFunction(26);
tr.SetFanoFactor(0.179);
// Connect the track to a sensor.
tr.SetSensor(&sensor);
DriftLineRKF drift;
drift.SetSensor(&sensor);
drift.SetMaximumStepSize(10.e-4);
// Plot the track and the drift lines.
ViewDrift driftView;
tr.EnablePlotting(&driftView);
drift.EnablePlotting(&driftView);
// Simulate an ion track.
tr.NewTrack(0., 0., 0., 0., 0., 1., 0.);
// Loop over the clusters.
double xc, yc, zc, tc, ec, ekin;
int ne = 0;
while (tr.GetCluster(xc, yc, zc, tc, ne, ec, ekin)) {
// Simulate electron and ion drift lines starting
// from the cluster position.
// Scale the induced current by the number of electron/ion pairs
// in the cluster.
drift.SetElectronSignalScalingFactor(ne);
drift.DriftElectron(xc, yc, zc, tc);
drift.SetHoleSignalScalingFactor(ne);
drift.DriftHole(xc, yc, zc, tc);
}
driftView.SetArea(-2.e-4, 0., 2.e-4, 100.e-4);
driftView.Plot(true);
ViewSignal signalView;
signalView.SetSensor(&sensor);
signalView.PlotSignal("readout", true, true, true);
app.Run();
//return 0;
}