Simulation of Fe55 source on a parallel plate detector

Dear Garfield experts,

I am using the fe55.C example modified in order to simulate the response of a micromegas detector on a point Fe55 iron X-ray source. Then I track the electrons from the photon absorption and calculate the deposited charge on the mesh using the RKF method at first.
In the left picture (number of electrons) are the electrons produced which seems to be in the right relative position concerning the Argon escape peak and the Iron dominant peak.
In the middle picture (number of electrons 2), where the number of electrons is shown on the mesh, why the argon (??) peak is placed too far from the dominant one?
The third picture on the right is the experimental one that we would like to simulate…

Any suggestion??

The code I used:

// X-ray conversion
// -------------------------------------------------------------------
// Simulate the Fe55 spectrum
//
// Thanks to Dorothea Pfeiffer and Heinrich Schindler from CERN for their help.
// -------------------------------------------------------------------
// Lucian Scharenberg
// scharenberg@physik.uni-bonn.de
// 05 APR 2018


#include <iostream>
#include <fstream>
#include <cmath>

#include <TCanvas.h>
#include <TROOT.h>
#include <TApplication.h>
#include <TH1F.h>

#include "Garfield/TrackHeed.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/SolidTube.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/FundamentalConstants.hh"
#include "Garfield/Random.hh"
#include "Garfield/Plotting.hh"
#include "Garfield/ComponentAnalyticField.hh"

#include "Garfield/ViewDrift.hh"
#include "Garfield/DriftLineRKF.hh"

#include "Garfield/AvalancheMicroscopic.hh"

using namespace Garfield;

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

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


  MediumMagboltz gas;  
  gas.LoadGasFile("Ar95_iso5_1atm.gas");
 
 
   // Build the Micromegas geometry.
  // We use separate components for drift and amplification regions. 

  // 50 um amplification gap. 
  ComponentAnalyticField cmpA; 
  cmpA.SetMedium(&gas);
  constexpr double yMesh = 0.005;
  constexpr double vMesh = 0.;
  constexpr double vPad = 350.;
  cmpA.AddPlaneY(yMesh, vMesh, "mesh1");
  cmpA.AddPlaneY(0., vPad, "pad");
  cmpA.AddReadout("mesh1");
  cmpA.AddReadout("pad");
  // 10 mm drift gap.
  ComponentAnalyticField cmpD;
  cmpD.SetMedium(&gas);
  constexpr double yDrift = yMesh + 1.0; 
  constexpr double vDrift = -600.; 
  cmpD.AddPlaneY(yDrift, vDrift, "drift");
  cmpD.AddPlaneY(yMesh, vMesh, "mesh2");
  cmpD.AddReadout("mesh2");
  cmpD.AddReadout("drift");
  // Assemble a sensor.
  Sensor sensor;
  sensor.AddComponent(&cmpA); 
  sensor.AddComponent(&cmpD); 
  sensor.SetArea(-5.0, 0., -5.0, 5.0, yDrift, 5.0);
  
  sensor.AddElectrode(&cmpA, "mesh1"); 
  sensor.AddElectrode(&cmpA, "pad"); 
  sensor.AddElectrode(&cmpD, "mesh2");
  sensor.AddElectrode(&cmpD, "drift");
 
 
   // RKF integration.
  DriftLineRKF drift;
  drift.SetSensor(&sensor);
 

  // Use Heed for simulating the photon absorption.
  TrackHeed track(&sensor);
  track.EnableElectricField();
  // Histogram
  const int nBins = 500;
  TH1::StatOverflows(true);
  TH1F hElectrons("hElectrons", ";Number of electrons;", 
                  nBins, -0.5, nBins - 0.5);
                  
  TH1F hElectrons2("hElectrons2", ";Number of electrons2;", 
                  200, -0.5, 80000 - 0.5);                 

                                   
  const unsigned int nEvents = 10000;

	  for (unsigned int i = 0; i < nEvents; ++i) {
   		 if (i % 1000 == 0) std::cout << i << "/" << nEvents << "\n";
    // Initial coordinates of the photon.

    int nsum = 0;
    int ne = 0;

    const double x0 = 0.;
    const double y0 = 1.0;
    const double z0 = 0.;
    const double t0 = 0.;

    const double egamma = 5900;

    track.TransportPhoton(x0, y0, z0, t0, egamma, 0., -1.0, 0., ne);

    	if (ne > 0) hElectrons.Fill(ne);

		for (int k = 0; k < ne; ++k) {

	
      nsum += ne;	


  double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.;
  double dx = 0., dy = 0., dz = 0.;
  track.GetElectron(k, xe, ye, ze, te, ee, dx, dy, dz);


  drift.DriftElectron(xe, ye, ze, te);


      double xee = 0., yee = 0., zee = 0., tee = 0.;
      int status = 0;
      drift.GetEndPoint(xee, yee, zee, tee, status);
      yee = yMesh - 1.e-10;
      drift.DriftElectron(xee, yee, zee, tee);

std::cout << "k=" << k+1 << " ne=" << ne << " xe=" << xe << " ye=" << ye << " ze=" << ze << " te=" << te << " xee=" << xee << " yee=" << yee << " zee=" << zee << " tee=" << tee << " nsum=" << nsum <<"\n";
}

    if (nsum > 0)     hElectrons2.Fill(nsum);

 	if (nsum > 0)      std::cout << "nsum: " << nsum << "\n";
 
  }

  TCanvas c("c", "", 600, 600);
  c.cd();
  hElectrons.SetFillColor(kBlue + 2);
  hElectrons.SetLineColor(kBlue + 2);
  hElectrons.Draw();
   c.Update(); 
   
   
   TCanvas c2("c2", "", 600, 600);
  c2.cd();
  hElectrons2.SetFillColor(kBlue + 2);
  hElectrons2.SetLineColor(kBlue + 2);
  hElectrons2.Draw(); 
  
  
  c2.Update();
  app.Run(true);
}

In your internal loop for (int k = 0; k < ne; ++k) you are adding ne to nsum ne times, which basically means you are getting square of ne. Instead you want to add drift.GetNumberOfElectronEndpoints() to nsum or calculate signals (Signals | Garfield++) and integrate. Also, not sure why you are drifting electrons twice.

Dear expert,

Thanks so much for your suggestion.
I commented the last drift.electron statement as you pointed out.
Unfortunately the **drift.GetNumberOfElectronEndpoints()** is not compatible with the RKF class I am using to drift the produced by the X-rays electrons to the mesh where I would like to collect them (at yee=0.005). So, I used instead the condition: nsum += drift.GetNumberOfDriftLinePoints();
The results are in the right direction, but I doubt if I am correct…
I wonder whether you could have a look and let me know.

Best regards

Georgios

// X-ray conversion
// -------------------------------------------------------------------
// Simulate the Fe55 spectrum
//
// Thanks to Dorothea Pfeiffer and Heinrich Schindler from CERN for their help.
// -------------------------------------------------------------------
// Lucian Scharenberg
// scharenberg@physik.uni-bonn.de
// 05 APR 2018

#include <iostream>
#include <fstream>
#include <cmath>

#include <TCanvas.h>
#include <TROOT.h>
#include <TApplication.h>
#include <TH1F.h>

#include "Garfield/TrackHeed.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/SolidTube.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/FundamentalConstants.hh"
#include "Garfield/Random.hh"
#include "Garfield/Plotting.hh"
#include "Garfield/ComponentAnalyticField.hh"

#include "Garfield/ViewDrift.hh"
#include "Garfield/DriftLineRKF.hh"

#include "Garfield/AvalancheMicroscopic.hh"

using namespace Garfield;

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

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

  MediumMagboltz gas;  
  gas.LoadGasFile("Ar95_iso5_1atm.gas");
 
 
   // Build the Micromegas geometry.
  // We use separate components for drift and amplification regions. 
  // 50 um amplification gap. 
  
  ComponentAnalyticField cmpA; 
  cmpA.SetMedium(&gas);
  constexpr double yMesh = 0.005;
  constexpr double vMesh = 0.;
  constexpr double vPad = 350.;
  cmpA.AddPlaneY(yMesh, vMesh, "mesh1");
  cmpA.AddPlaneY(0., vPad, "pad");
  cmpA.AddReadout("mesh1");
  cmpA.AddReadout("pad");
  // 10 mm drift gap.
  ComponentAnalyticField cmpD;
  cmpD.SetMedium(&gas);
  constexpr double yDrift = yMesh + 1.0; 
  constexpr double vDrift = -600.; 
  cmpD.AddPlaneY(yDrift, vDrift, "drift");
  cmpD.AddPlaneY(yMesh, vMesh, "mesh2");
  cmpD.AddReadout("mesh2");
  cmpD.AddReadout("drift");
  // Assemble a sensor.
  Sensor sensor;
  sensor.AddComponent(&cmpA); 
  sensor.AddComponent(&cmpD); 
  sensor.SetArea(-5.0, 0., -5.0, 5.0, yDrift, 5.0);
  
  sensor.AddElectrode(&cmpA, "mesh1"); 
  sensor.AddElectrode(&cmpA, "pad"); 
  sensor.AddElectrode(&cmpD, "mesh2");
  sensor.AddElectrode(&cmpD, "drift");
 
 
   // RKF integration.
  DriftLineRKF drift;
  drift.SetSensor(&sensor);
 
//  AvalancheMicroscopic aval;  
//  aval.SetSensor(&sensor);    

  // Use Heed for simulating the photon absorption.
  TrackHeed track(&sensor);
  track.EnableElectricField();
  // Histogram
  const int nBins = 500;
  TH1::StatOverflows(true);
  TH1F hElectrons("hElectrons", ";Number of electrons;", 
                  nBins, -0.5, nBins - 0.5);
                  
   TH1F hElectrons2("hElectrons2", ";Number of electrons2;", 
                  200, -0.5, 10000 - 0.5);                 

                                   
  const unsigned int nEvents = 10000;
  for (unsigned int i = 0; i < nEvents; ++i) {
    if (i % 1000 == 0) std::cout << i << "/" << nEvents << "\n";
  

    int nsum = 0;
    int ne = 0;

    // Initial coordinates of the photon.
    const double x0 = 0.;
    const double y0 = 1.0;
    const double z0 = 0.;
    const double t0 = 0.;
//    Sample the photon energy, using the relative intensities according to XDB.
//    const double r = 167. * RndmUniform();
//    const double egamma = r < 100. ? 5898.8 : r < 150. ? 5887.6 : 6490.4; 
    const double egamma = 5900;  

    track.TransportPhoton(x0, y0, z0, t0, egamma, 0., -1.0, 0., ne);

    if (ne > 0) hElectrons.Fill(ne);


	for (int k = 0; k < ne; ++k) {
     
      nsum += drift.GetNumberOfDriftLinePoints();
    

  double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.;
  double dx = 0., dy = 0., dz = 0.;

  track.GetElectron(k, xe, ye, ze, te, ee, dx, dy, dz);
  drift.DriftElectron(xe, ye, ze, te);

      double xee = 0., yee = 0., zee = 0., tee = 0.;
      int status = 0;
      drift.GetEndPoint(xee, yee, zee, tee, status);
      yee = yMesh - 1.e-10;
//      drift.DriftElectron(xee, yee, zee, tee);

std::cout << " k=" << k+1 << " ne=" << ne << " DriftLInePoints=" << drift.GetNumberOfDriftLinePoints() << " xe=" << xe << " ye=" << ye << " ze=" << ze << " te=" << te << " xee=" << xee << " yee=" << yee << " zee=" << zee << " tee=" << tee << " nsum=" << nsum <<"\n";
}

    if (nsum > 0)     hElectrons2.Fill(nsum);
 
    if (nsum > 0)      std::cout << "nsum: " << nsum << "\n";
 
  }

  TCanvas c("c", "", 600, 600);
  c.cd();
  hElectrons.SetFillColor(kBlue + 2);
  hElectrons.SetLineColor(kBlue + 2);
  hElectrons.Draw();
   c.Update(); 
   
   
   TCanvas c2("c2", "", 600, 600);
  c2.cd();
  hElectrons2.SetFillColor(kBlue + 2);
  hElectrons2.SetLineColor(kBlue + 2);
  hElectrons2.Draw(); 
  
  
  c2.Update();
  app.Run(true);
}

Ah, it’s RKF, right. DriftLineRKF does not simulate avalanches, just a single drift line. GetNumberOfDriftLinePoints() will give you number of points (steps) in that line, so that’s wrong. DriftLineRKF can however obtain the number of electrons at the end of would-be avalanche by integrating Townsend formula - you can extract that using GetAvalancheSize() method.

// X-ray conversion
// -------------------------------------------------------------------
// Simulate the Fe55 spectrum
//
// Thanks to Dorothea Pfeiffer and Heinrich Schindler from CERN for their help.
// -------------------------------------------------------------------
// Lucian Scharenberg
// scharenberg@physik.uni-bonn.de
// 05 APR 2018

#include <iostream>
#include <fstream>
#include <cmath>

#include <TCanvas.h>
#include <TROOT.h>
#include <TApplication.h>
#include <TH1F.h>

#include "Garfield/TrackHeed.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/SolidTube.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/FundamentalConstants.hh"
#include "Garfield/Random.hh"
#include "Garfield/Plotting.hh"
#include "Garfield/ComponentAnalyticField.hh"

#include "Garfield/ViewDrift.hh"
#include "Garfield/DriftLineRKF.hh"

#include "Garfield/AvalancheMicroscopic.hh"

using namespace Garfield;

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

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

  MediumMagboltz gas;  
  gas.LoadGasFile("Ar95_iso5_1atm.gas");
 
 
   // Build the Micromegas geometry.
  // We use separate components for drift and amplification regions. 
  // 50 um amplification gap. 
  
  ComponentAnalyticField cmpA; 
  cmpA.SetMedium(&gas);
  constexpr double yMesh = 0.005;
  constexpr double vMesh = 0.;
  constexpr double vPad = 350.;
  cmpA.AddPlaneY(yMesh, vMesh, "mesh1");
  cmpA.AddPlaneY(0., vPad, "pad");
  cmpA.AddReadout("mesh1");
  cmpA.AddReadout("pad");
  // 10 mm drift gap.
  ComponentAnalyticField cmpD;
  cmpD.SetMedium(&gas);
  constexpr double yDrift = yMesh + 1.0; 
  constexpr double vDrift = -600.; 
  cmpD.AddPlaneY(yDrift, vDrift, "drift");
  cmpD.AddPlaneY(yMesh, vMesh, "mesh2");
  cmpD.AddReadout("mesh2");
  cmpD.AddReadout("drift");
  // Assemble a sensor.
  Sensor sensor;
  sensor.AddComponent(&cmpA); 
  sensor.AddComponent(&cmpD); 
  sensor.SetArea(-5.0, 0., -5.0, 5.0, yDrift, 5.0);
  
  sensor.AddElectrode(&cmpA, "mesh1"); 
  sensor.AddElectrode(&cmpA, "pad"); 
  sensor.AddElectrode(&cmpD, "mesh2");
  sensor.AddElectrode(&cmpD, "drift");
 
 
   // RKF integration.
//  DriftLineRKF drift;
//  drift.SetSensor(&sensor);
 
   AvalancheMicroscopic drift;  
   drift.SetSensor(&sensor);	

  // Use Heed for simulating the photon absorption.
  TrackHeed track(&sensor);
  track.EnableElectricField();
  // Histogram
  const int nBins = 500;
  TH1::StatOverflows(true);
  TH1F hElectrons("hElectrons", ";Number of electrons;", 
                  200, -0.5, nBins - 0.5);
                  
   TH1F hElectrons2("hElectrons2", ";Number of electrons2;", 
                  200, -0.5, nBins - 0.5);                 

                                   
  const unsigned int nEvents = 100;
  for (unsigned int i = 0; i < nEvents; ++i) {
    if (i % 10 == 0) std::cout << i << "/" << nEvents << "\n";
  

    int nsum = 0;
    int ne = 0;

    // Initial coordinates of the photon.
    const double x0 = 0.;
    const double y0 = 1.0;
    const double z0 = 0.;
    const double t0 = 0.;
//    Sample the photon energy, using the relative intensities according to XDB.
    const double r = 167. * RndmUniform();
    const double egamma = r < 100. ? 5898.8 : r < 150. ? 5887.6 : 6490.4; 
//    const double egamma = 5900;  

    track.TransportPhoton(x0, y0, z0, t0, egamma, 0., -1.0, 0., ne);

    if (ne > 0) hElectrons.Fill(ne);


	for (int k = 0; k < ne; ++k) {
     
      nsum += drift.GetNumberOfElectronEndpoints();
    

  double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.1;
  double dx = 0., dy = 0., dz = 0.;

  track.GetElectron(k, xe, ye, ze, te, ee, dx, dy, dz);
//  drift.DriftElectron(xe, ye, ze, te);
          drift.DriftElectron(xe, ye, ze, te, 0.1, 0., 0., 0.);

          // Get the electron's end point.
          double x1 = 0., y1 = 0., z1 = 0., t1 = 0., e1 = 0.;
          int status = 0;
  ye = yMesh - 1.e-10;
          drift.GetElectronEndpoint(0, x1, y1, z1, t1, e1, 
                                      xe, ye, ze, te, ee, status);
//      drift.GetEndPoint(xee, yee, zee, tee, status);
//      yee = yMesh - 1.e-10;
//    drift.GetElectronEndpoint(0, xee, yee, z1, t1, e1, xe, ye, ze, te, ee, status);
//      drift.DriftElectron(xee, yee, zee, tee);


std::cout << " k=" << k+1 << " ne=" << ne << " NumberOfElectronEndpoints=" << drift.GetNumberOfElectronEndpoints() << " x1=" << x1 << " y1=" << y1 << " z1=" << z1 << " t1=" << t1 << " xe=" << xe << " ye=" << ye << " ze=" << ze << " te=" << te << " nsum=" << nsum <<"\n";


//std::cout << " k=" << k+1 << " ne=" << ne << " NumberOfElectronEndpoints=" << drift.GetNumberOfElectronEndpoints() << " xe=" << xe << " ye=" << ye << " ze=" << ze << " te=" << te << " xee=" << xee << " yee=" << yee << " zee=" << zee << " tee=" << tee << " nsum=" << nsum <<"\n";
}

    if (nsum > 0)     hElectrons2.Fill(nsum);
 
 //   if (nsum > 0)      std::cout << "nsum: " << nsum << "\n";
 
  }

  TCanvas c("c", "", 600, 600);
  c.cd();
  hElectrons.SetFillColor(kBlue + 2);
  hElectrons.SetLineColor(kBlue + 2);
  hElectrons.Draw();
   c.Update(); 
   
   
   TCanvas c2("c2", "", 600, 600);
  c2.cd();
  hElectrons2.SetFillColor(kBlue + 2);
  hElectrons2.SetLineColor(kBlue + 2);
  hElectrons2.Draw(); 
  
  
  c2.Update();
  app.Run(true);
}

Dear sir,

I have just tried to work on the code attached using the AvalancheMicroscopic method.
As you could see in the picts, with just 100 events (much slower the RKF as expected) both plots are identical (??) and the charge deposited on the escape peak is not spread alot… Is it just a matter of statistics or more mistakes are present?

Best regards

Georgios

This looks fine for low statistics, it does not look like there are other mistakes. Try to run with decent statistics (5k or 10k events).
Kind regards
Piet

Hi,
AvalancheMicroscopic::DriftElectron simulates the drift line of a single electron without tracking the secondary electrons in the avalanche (so GetNumberOfElectronEndpoints should always return 1). If you want to simulate an avalanche, you should use AvalancheMicroscopic::AvalancheElectron. By the way, you need to call GetNumberOfElectronEndpoints after simulating the electron drift or avalanche, not before.

track.TransportPhoton(x0, y0, z0, t0, egamma, 0., -1.0, 0., ne);
if (ne > 0) hElectrons.Fill(ne);
for (int k = 0; k < ne; ++k) {
  double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.1;
  double dx = 0., dy = 0., dz = 0.;
  track.GetElectron(k, xe, ye, ze, te, ee, dx, dy, dz);
  drift.AvalancheElectron(xe, ye, ze, te, 0.1, 0., 0., 0.);
  // Get the number of electrons in the avalanche.
  nsum += drift.GetElectrons().size();
}

But you can also use DriftLineRKF and retrieve the avalanche size as @mjag suggested.


Hi, thanks so much for your remarks. I tried to follow your suggestion but in vain…
It seems that both plots are identical and also the drift.GetElectrons().size() equals 1.
I wonder where is the mistake… I was expecting that the right plot could reproduce the experimental one.

Best regards

// X-ray conversion
// -------------------------------------------------------------------
// Simulate the Fe55 spectrum
//
// Thanks to Dorothea Pfeiffer and Heinrich Schindler from CERN for their help.
// -------------------------------------------------------------------
// Lucian Scharenberg
// scharenberg@physik.uni-bonn.de
// 05 APR 2018

#include <iostream>
#include <fstream>
#include <cmath>

#include <TCanvas.h>
#include <TROOT.h>
#include <TApplication.h>
#include <TH1F.h>

#include "Garfield/TrackHeed.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/SolidTube.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/FundamentalConstants.hh"
#include "Garfield/Random.hh"
#include "Garfield/Plotting.hh"
#include "Garfield/ComponentAnalyticField.hh"

#include "Garfield/ViewDrift.hh"
#include "Garfield/DriftLineRKF.hh"

#include "Garfield/AvalancheMicroscopic.hh"

using namespace Garfield;

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

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

 MediumMagboltz gas;  
 gas.LoadGasFile("Ar95_iso5_1atm.gas");


  // Build the Micromegas geometry.
 // We use separate components for drift and amplification regions. 
 // 50 um amplification gap. 
 
 ComponentAnalyticField cmpA; 
 cmpA.SetMedium(&gas);
 constexpr double yMesh = 0.005;
 constexpr double vMesh = 0.;
 constexpr double vPad = 350.;
 cmpA.AddPlaneY(yMesh, vMesh, "mesh1");
 cmpA.AddPlaneY(0., vPad, "pad");
 cmpA.AddReadout("mesh1");
 cmpA.AddReadout("pad");
 // 10 mm drift gap.
 ComponentAnalyticField cmpD;
 cmpD.SetMedium(&gas);
 constexpr double yDrift = yMesh + 1.0; 
 constexpr double vDrift = -600.; 
 cmpD.AddPlaneY(yDrift, vDrift, "drift");
 cmpD.AddPlaneY(yMesh, vMesh, "mesh2");
 cmpD.AddReadout("mesh2");
 cmpD.AddReadout("drift");
 // Assemble a sensor.
 Sensor sensor;
 sensor.AddComponent(&cmpA); 
 sensor.AddComponent(&cmpD); 
 //sensor.SetArea(-5.0, 0., -5.0, 5.0, yDrift, 5.0);
 sensor.SetArea(-3.0, 0., -3.0, 3.0, yDrift, 3.0);  
 sensor.AddElectrode(&cmpA, "mesh1"); 
 sensor.AddElectrode(&cmpA, "pad"); 
 sensor.AddElectrode(&cmpD, "mesh2");
 sensor.AddElectrode(&cmpD, "drift");


  // RKF integration.
//  DriftLineRKF drift;
//  drift.SetSensor(&sensor);

  AvalancheMicroscopic drift;  
  drift.SetSensor(&sensor);	

 // Use Heed for simulating the photon absorption.
 TrackHeed track(&sensor);
 track.EnableElectricField();
 // Histogram
 const int nBins = 500;
 TH1::StatOverflows(true);
 TH1F hElectrons("hElectrons", ";Number of electrons;", 
                 200, -0.5, nBins - 0.5);
                 
  TH1F hElectrons2("hElectrons2", ";Number of electrons2;", 
                 200, -0.5, nBins - 0.5);                 

                                  
 const unsigned int nEvents = 100;
 for (unsigned int i = 0; i < nEvents; ++i) {
   if (i % 10 == 0) std::cout << i << "/" << nEvents << "\n";
 

   int nsum = 0;
   int ne = 0;

   // Initial coordinates of the photon.
   const double x0 = 0.;
   const double y0 = 1.0;
   const double z0 = 0.;
   const double t0 = 0.;
//    Sample the photon energy, using the relative intensities according to XDB.
   const double r = 167. * RndmUniform();
   const double egamma = r < 100. ? 5898.8 : r < 150. ? 5887.6 : 6490.4; 
//    const double egamma = 5900;  

   track.TransportPhoton(x0, y0, z0, t0, egamma, 0., -1.0, 0., ne);

   if (ne > 0) hElectrons.Fill(ne);

   for (int k = 0; k < ne; ++k) {
    
//      nsum += drift.GetNumberOfElectronEndpoints();
   
 double xe = 0., ye = 0., ze = 0., te = 0., ee = 0.1;
 double dx = 0., dy = 0., dz = 0.;

 track.GetElectron(k, xe, ye, ze, te, ee, dx, dy, dz);

//  drift.DriftElectron(xe, ye, ze, te);
//          drift.DriftElectron(xe, ye, ze, te, 0.1, 0., 0., 0.);

 drift.AvalancheElectron(xe, ye, ze, te, 0.1, 0., 0., 0.);
 nsum += drift.GetElectrons().size();
         // Get the electron's end point.
         double x1 = 0., y1 = 0., z1 = 0., t1 = 0., e1 = 0.;
         int status = 0;
 ye = yMesh - 1.e-10;
         drift.GetElectronEndpoint(0, x1, y1, z1, t1, e1, 
                                     xe, ye, ze, te, ee, status);
//      drift.GetEndPoint(xee, yee, zee, tee, status);
//      yee = yMesh - 1.e-10;
//    drift.GetElectronEndpoint(0, xee, yee, z1, t1, e1, xe, ye, ze, te, ee, status);
//      drift.DriftElectron(xee, yee, zee, tee);


std::cout << " i=" << i << " k=" << k+1 << " ne=" << ne << " NumberOfElectronEndpoints=" << drift.GetNumberOfElectronEndpoints() <<" size=" << drift.GetElectrons().size() << " x1=" << x1 << " y1=" << y1 << " z1=" << z1 << " t1=" << t1 << " xe=" << xe << " ye=" << ye << " ze=" << ze << " te=" << te << " nsum=" << nsum <<"\n";


//std::cout << " k=" << k+1 << " ne=" << ne << " NumberOfElectronEndpoints=" << drift.GetNumberOfElectronEndpoints() << " xe=" << xe << " ye=" << ye << " ze=" << ze << " te=" << te << " xee=" << xee << " yee=" << yee << " zee=" << zee << " tee=" << tee << " nsum=" << nsum <<"\n";
}

   if (nsum > 0)     hElectrons2.Fill(nsum);

//   if (nsum > 0)      std::cout << "nsum: " << nsum << "\n";

 }

 TCanvas c("c", "", 600, 600);
 c.cd();
 hElectrons.SetFillColor(kBlue + 2);
 hElectrons.SetLineColor(kBlue + 2);
 hElectrons.Draw();
  c.Update(); 
  
  
  TCanvas c2("c2", "", 600, 600);
 c2.cd();
 hElectrons2.SetFillColor(kBlue + 2);
 hElectrons2.SetLineColor(kBlue + 2);
 hElectrons2.Draw(); 
 
 
 c2.Update();
 app.Run(true);
}

OK, I see why electron was drifted twice in original code. The code simplifies the mesh with an infinite plate. Now, since it’s a solid plate, electrons stop on it, so they need to be moved slightly below the “mesh” and simulated it again.

To generalize it for AvalancheMicroscopic (which I assume is what you want to do) adds some complexity, since, in rare cases, the electrons can actually multiply in drift (especially if you increase drift voltage). Basically you need to loop over the electrons that arrive at the mesh and resimulate all of them:

#include <iostream>
#include <fstream>
#include <cmath>

#include <TCanvas.h>
#include <TROOT.h>
#include <TApplication.h>
#include <TH1F.h>

#include "Garfield/TrackHeed.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/SolidTube.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/ComponentConstant.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/FundamentalConstants.hh"
#include "Garfield/Random.hh"
#include "Garfield/Plotting.hh"
#include "Garfield/ComponentAnalyticField.hh"

#include "Garfield/ViewDrift.hh"
#include "Garfield/DriftLineRKF.hh"

#include "Garfield/AvalancheMicroscopic.hh"

using namespace Garfield;

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

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

    MediumMagboltz gas;
    gas.LoadGasFile("Ar95_iso5_1atm.gas");

    // Build the Micromegas geometry.
    // We use separate components for drift and amplification regions.
    // 50 um amplification gap.
    ComponentAnalyticField cmpA;
    cmpA.SetMedium(&gas);
    constexpr double yMesh = 0.005;
    constexpr double vMesh = 0.;
    constexpr double vPad = 350.;
    cmpA.AddPlaneY(yMesh, vMesh, "mesh1");
    cmpA.AddPlaneY(0., vPad, "pad");
    cmpA.AddReadout("mesh1");
    cmpA.AddReadout("pad");
    Sensor sensorA;
    sensorA.AddComponent(&cmpA);
    sensorA.SetArea(-3.0, 0., -3.0, 3.0, yMesh, 3.0);

    // 10 mm drift gap.
    ComponentAnalyticField cmpD;
    cmpD.SetMedium(&gas);
    constexpr double yDrift = yMesh + 1.0;
    constexpr double vDrift = -600.;
    cmpD.AddPlaneY(yDrift, vDrift, "drift");
    cmpD.AddPlaneY(yMesh, vMesh, "mesh2");
    cmpD.AddReadout("mesh2");
    cmpD.AddReadout("drift");
    Sensor sensorD;
    sensorD.AddComponent(&cmpD);
    sensorD.SetArea(-3.0, yMesh, -3.0, 3.0, yDrift, 3.0);

    AvalancheMicroscopic driftA, driftD;
    driftA.SetSensor(&sensorA);
    driftD.SetSensor(&sensorD);

    // Use Heed for simulating the photon absorption.
    TrackHeed track(&sensorD);
    track.EnableElectricField();

    // Histogram
    const int nBins = 500;
    TH1::StatOverflows(true);
    TH1F hElectrons("hElectrons", ";Number of electrons;",
                    nBins, -0.5, nBins - 0.5);

    TH1F hElectrons2("hElectrons2", ";Number of electrons2;",
                    200, -0.5, 80000 - 0.5);

    // Simulate 100 events
    const unsigned int nEvents = 100;
    for (unsigned int i = 0; i < nEvents; ++i) {
        if (i % 10 == 0) std::cout << i << "/" << nEvents << "\n";
        int nsum = 0;
        int ne = 0;

        // Initial coordinates of the photon.
        const double x0 = 0.;
        const double y0 = 1.0;
        const double z0 = 0.;
        const double t0 = 0.;
        // Sample the photon energy, using the relative intensities according to XDB.
        const double r = 167. * RndmUniform();
        const double egamma = r < 100. ? 5898.8 : r < 150. ? 5887.6 : 6490.4;

        track.TransportPhoton(x0, y0, z0, t0, egamma, 0., -1.0, 0., ne);

        // Fill number of electrons in ionization
        if (ne > 0) hElectrons.Fill(ne);

        // Loop over electrons in ionization
        for (int k1 = 0; k1 < ne; ++k1) {
            double xe, ye, ze, te, ee, dx, dy, dz;
            track.GetElectron(k1, xe, ye, ze, te, ee, dx, dy, dz);

            // Avalanche electrons in drifting area
            // (should not create much of an avalanche)
            driftD.AvalancheElectron(xe, ye, ze, te, 0, 0, 0, 0);

            // Loop over electrons that arrived at the mesh
            int n_mesh = driftD.GetNumberOfElectronEndpoints();
            for (int k2 = 0; k2 < n_mesh; ++k2) {
                double x1, y1, z1, t1, e1, x2, y2, z2, e2, t2;
                int status;
                driftD.GetElectronEndpoint(k2, x1, y1, z1, t1, e1, x2, y2, z2, t2, e2, status);

                // Resimulate - actual avalanche
                driftA.AvalancheElectron(x2, yMesh - 1e-4, z2, t2, 0, 0, 0, 0);

                int n_pad = driftA.GetNumberOfElectronEndpoints();
                nsum += n_pad;
            }
        }
        // Fill number of electrons at the pad
        if (nsum > 0) hElectrons2.Fill(nsum);
    }

    TCanvas c("c", "", 600, 600);
    c.cd();
    hElectrons.SetFillColor(kBlue + 2);
    hElectrons.SetLineColor(kBlue + 2);
    hElectrons.Draw();
    c.Update();

    TCanvas c2("c2", "", 600, 600);
    c2.cd();
    hElectrons2.SetFillColor(kBlue + 2);
    hElectrons2.SetLineColor(kBlue + 2);
    hElectrons2.Draw();

    c2.Update();
    app.Run(true);

Also see Examples/Asacusa/example.cpp · master · garfield / garfieldpp · GitLab, which does basically the same, except with some more setup, pions and assumes single track in drift.

Of course, since it’s a simplified geometry, it’s all unnecessary. You will note with microscopic method it is quite slow. You can use the solution with RKF using GetAvalancheSize() like I mentioned in my second reply. Microscopic would become necessary if you wanted to simulate small scale interactions with mesh.

Dear expert,

Thanks so much for your post. I got your code, and run it for 6000 events that lasted ~ 10h!
It is really quite slow as you pointed out but the right plot represents quite well the experimental amplitude distribution…


The results show the #electrons as given by the x-axis taken from the Mean is ~215 initially and ~116000 in the pad, leading to a ratio ~ 540. 215 e- in the ionization sounds goods (i.e about 27 eV needed per electron created) but I was expecting a ratio about 10 times bigger. 116000 e- corresponds to ~18.5 fC while I think one measured about 200 fC!.
I was thinking that might the driftA.AvalancheElectron(x2, yMesh - 1e-4, z2, t2, 0, 0, 0, 0) should be substituted with driftA.AvalancheElectron(x2, yMesh - 1.e-6, z2, t2, 0, 0, 0, 0) but this leads to much longer execution time and the ratio is improved by a little (663 vs 540)…

I tried to think about using the RKF as you suggested using aval.GetAvalancheSize(ne, ni) but I wasnt able to find any example using this solution with a method different than the AvalancheMicroscopic… I wonder whether you could give me some help…

Best regards

Dear @gtsileda

You can find an example using DriftLineRKF here. All available methods of the DriftLinkeRKF class are located here.

In short you should
(1) initialize an instance the class and link the sensor to this object:

DriftLineRKF drift;
drift.SetSensor(&sensor);

(2) drift all your electrons:

drift.DriftElectron(electron.x, electron.y, electron.z, electron.t);

(3) obtain the gain:

drift.GetAvalancheSize(ne, ni);

Good luck and let us know if you succeeded!
greets
Piet

Hello again. That shift should be as small as possible without causing numerical edge case problems. I changed it to 1e-4 for a test and forgot to change it back to 1e-10 as in original code, my bad.

As for RKF, simply adapt your original code:

  1. Remove the line with nsum += ne;
  2. Get avalanche size after drift.DriftElectron(xee, yee, zee, tee);:
double ne_pad, ni_mesh;
drift.GetAvalancheSize(ne_pad, ni_mesh);
nsum += ne_pad;
  1. Change upper limit for hElectrons2 to something appropriate for your gas.

As Piet suggested, Doxygen is a good source to find reference for different methods.

Best,
Michał