Help: Garfield++ :: AvalancheMicroscopic and TrackHeed

Hello everyone,

I am working on a simulation with Garfield++ where I want to analyze the microscopic details of the ionization process, particularly focusing on the secondary ionization generated during the avalanche process. I want to ensure that all electrons, including those produced by ionization (delta electrons), are correctly tracked and contribute to the signal.

(Picture: I am attempting to simulate the drift and avalanche of electrons produced along a track. While the drift lines are visualized correctly, I am having issues with the secondary ionization process and how to properly collect and visualize avalanche electrons.)

Here’s a snippet of my code:

#include <iostream>
#include <cmath>
#include <TApplication.h>
#include <TCanvas.h>
#include "Garfield/ViewField.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/Plotting.hh"

using namespace Garfield;

int main(int argc, char * argv[])
{

  TApplication app("app", &argc, argv);
  plottingEngine.SetDefaultStyle();

  // Set the gas mixture.
  MediumMagboltz gas;
  gas.SetComposition("ar", 93., "co2", 7.);
//    gas.SetComposition("ar", 90., "c4h10", 10.);

  // Build the Micromegas geometry.
  // We use separate components for drift and amplification regions. 
  ComponentAnalyticField cmpD;
  ComponentAnalyticField cmpA; 
  cmpD.SetMedium(&gas);
  cmpA.SetMedium(&gas);
  // Set the magnetic field [Tesla].
  constexpr double bfield = 0.;
  cmpD.SetMagneticField(0, 0, bfield);
  cmpA.SetMagneticField(0, 0, bfield);
  

  // 100 um amplification gap. 
  constexpr double yMesh = 0.01;
  constexpr double vMesh = -500.0;
  cmpA.AddPlaneY(yMesh, vMesh, "mesh");
  cmpA.AddPlaneY(0.0, 0.0, "bottom");

  // 3 mm drift gap.
  constexpr double yDrift = yMesh + 0.3;   // cm
  constexpr double vDrift = -800.0;   // Volt
  cmpD.AddPlaneY(yDrift, vDrift, "top");
  cmpD.AddPlaneY(yMesh, vMesh, "mesh");

  // We place three strips along z direction on the anode plane
  // in order to be able to resolve the avalanche signal.
  // Strip pitch [cm].
  constexpr double pitch = 0.05;
  // Distance between two strips [cm].
  constexpr double interstrip = 0.01; 
  // Half-width of a strip.
  constexpr double hw = 0.5 * (pitch - interstrip);
  const double xStrip1 = hw;
  cmpA.AddStripOnPlaneY('z', 0., xStrip1 - hw, xStrip1 + hw, "strip1");
  const double xStrip2 = xStrip1 + pitch;
  cmpA.AddStripOnPlaneY('z', 0., xStrip2 - hw, xStrip2 + hw, "strip2");
  const double xStrip3 = xStrip2 + pitch;
  cmpA.AddStripOnPlaneY('z', 0., xStrip3 - hw, xStrip3 + hw, "strip3");
  const double xStrip4 = xStrip3 + pitch;
  cmpA.AddStripOnPlaneY('z', 0., xStrip4 - hw, xStrip4 + hw, "strip4");

  // We want to calculate the signals induced on the strips. 
  cmpA.AddReadout("strip1");
  cmpA.AddReadout("strip2");
  cmpA.AddReadout("strip3");
  cmpA.AddReadout("strip4");

  // Assemble a sensor.
  Sensor sensor;
  sensor.AddComponent(&cmpD); 
  sensor.AddComponent(&cmpA); 
  sensor.SetArea(0., 0., -1., 6 * pitch, yDrift, 1.);

  // Request signal calculation for the strip electrodes.
  sensor.AddElectrode(&cmpA, "strip1"); 
  sensor.AddElectrode(&cmpA, "strip2"); 
  sensor.AddElectrode(&cmpA, "strip3");
  sensor.AddElectrode(&cmpA, "strip4");

  // Set the time window [ns] for the signal calculation.
  const double tMin = 0.; 
  const double tMax = 100.; 
  const double tStep = 0.05; 
  const int nTimeBins = int((tMax - tMin) / tStep); 
  sensor.SetTimeWindow(0., tStep, nTimeBins);

  // We use microscopic tracking for simulating the electron avalanche.
  AvalancheMicroscopic aval;
  aval.SetSensor(&sensor); 
  // Switch on signal calculation. 
  //aval.EnableSignalCalculation(); 
  aval.EnableSignalCalculation(true); 
  aval.UseWeightingPotential(true);
  aval.EnableWeightingFieldIntegration(true);
  aval.UseInducedCharge(true);
  aval.EnableDriftLines(true);
  aval.EnableMagneticField();
  aval.SetElectronTransportCut(0.0); 
  
  // Simulate an ionizing particle (negative pion) using Heed.
  TrackHeed track;
  track.SetParticle("mu-");
  constexpr double momentum =300.e6; // [eV / c]
  track.SetMomentum(momentum);
  track.SetSensor(&sensor);
  track.EnableMagneticField();
  // track.EnableElectricField();
  track.EnableDeltaElectronTransport(); 
  track.Initialise(&gas, true);


  // Construct a viewer to visualise the drift lines.
  ViewDrift driftView;
  track.EnablePlotting(&driftView);
  aval.EnablePlotting(&driftView);
  aval.EnableExcitationMarkers(true); 
  aval.EnableIonisationMarkers(true); 
  aval.EnableAttachmentMarkers(true);


 // Set the starting point of the incoming ionising track.
  double xt = 0.1; // [cm]
  double yt = yDrift;
  double zt = 0.0;

  // Now simulate a track.
  track.NewTrack(xt, yt, zt, 0.87, 0.5, -1, 0);
  // Loop over the clusters.
  double xc, yc, zc, tc, ec, extra;
  int nc;
  while (track.GetCluster(xc, yc, zc, tc, nc, ec, extra)) 
  {
    for (int j = 0; j < nc; ++j) 
    {
      double xe, ye, ze, te, ee, dxe, dye, dze;
      track.GetElectron(j, xe, ye, ze, te, ee, dxe, dye, dze);
      // Simulate the drift/avalanche of this electron.
      aval.AvalancheElectron(xe, ye, ze, te, 100.0, dxe, dye, dze);
    }
  }

Hi,

Adding @hschindl in the loop.

D

1 Like

Dear Christina,

As soon as your electrons and ions encounter a plane during their trajectory they will be stoped. This means you will not be able to build up the toy-model of your micromegas in this way. In principle you just want to ‘glue’ two constant fields together, I would suggest doing it the following way:

 // Make a component for the drift field.
    ComponentUser cmp;
  
    // Parameterisation of the drift and amplification field.
    auto efield = [](const double /*x*/, const double y, const double /*z*/,
                     double& ex, double& ey, double& ez) {
        ex = ez = 0.;
        // Starting point of the gain layer.
        constexpr double y0 = 128.e-4;
        
        if (y <= y0) {
            ey = 40.e3; // V/cm
        } else {
            ey = 600.; // V/cm
        }
    };
    cmp.SetElectricField(efield);