The questions about the induced current

Garfield++
1

I want to get the induced current signal, however the result show in the below, I don’t know what’s the reason for the error and I modify the code several times, it is still the same.

image
image1607×850 108 KB

the code is in the below

#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 = 1000;
  const double tmin =  0.;
  const double tmax = 100.;
  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);

  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.1, 0.204, -0.07, 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(-1.67e-1, 0.17, -0.67e-1, 0.2);
  driftView.Plot(true);

  ViewSignal signalView;
  signalView.SetSensor(&sensor);
  signalView.PlotSignal("readout", true, true, true);

  app.Run();
  return 0;
}
2

As you can see from the output of the program, the import of the weighting field map failed:

ComponentAnsys123::SetWeightingField:
    Could not open file weight1.lis for reading.
    The file perhaps does not exist.

Make sure you copy the file to the directory from which you are running the executable.

3

Hi，
Thanks for your reply
I remember I have copied the weight1.lis file to the example directory’s subdirectory Gem.
I will check it again, and it must be a careless mistake.
Thanks for your help

4

Hi,
I have reproduce my programme, but the result is still the same as below

image
image1644×839 115 KB

would you please help me to find out the reason, Thanks a lot.

5

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).

6

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.

7

Hi,
I have deleted all the codes about the drift lines
the result is still the same in the below

image
image1301×755 103 KB

8

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;
}
9

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
10

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.

image
image1607×850 108 KB

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.

11

Besides, how to get the number of electrons/ions arriving at each electrode. Which code or output files can make it

12
  • 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 (.lis files).
  • 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 Initialise function (it’s the last argument of the function, in the code you posted earlier it was "cm").
13

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.

14

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;
  // ...
}
15

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.

16

Ok, so maybe I should set the length unit to mm

17

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);
}
18

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?

19

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;
}
20

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);
      }
    }
  }