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.

1 Like

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

1 Like

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ł

1 Like

Dear expert,

Thanks so much for your support. The RKF works following your suggestions and evaluate much faster the magnitude of the gain after avalanche. But I am afraid the only way to go is the much slower microscopic MC method where the diffusion takes place resulting in a spectra at the end that matches the experimental amplitude distribution. With RKF, both distributions (electrons after ionization and electrons at the pad) are statistically the same with the exception that the 2nd one has in the x-axis the number of electrons x gain - since there is no diffusion that could widen the 1st spectra.

Best regards

Georgios

Dear Georgios

If you would like to have the best possible agreement with data, you should indeed go with Microscopic Tracking, but you would also need to simulate the electric field in a 3D structure with your exact Mesh geometry (mesh wire thickness, lines/inch, mesh fabrication technique: woven/calendered/…) using for instance Finite Element Methods and import the solution in garfield. This would allow to estimate correctly the transversal spread and losses of electrons attached to the Mesh.

Simpler than Microscopic Tracking, but more advanced than DriftLineRFK would be the use of the AvalancheMC class that incorporates the longitudinal and transversal diffusion. Maybe this could be enough to have agreement with your data? Would you mind showing how good your data/simulation agreement is right now for DriftLineRKF and AvalancheMicroscopic?

greets
Piet

1 Like

Dear Piet,

Thanks a lot for your mail and your suggestions. Below are the plots for electrons after ionization and electrons collected in the pad with RKF, Microscopic MC methods and the experimental data. The resolution achieved with the latter is ~ 14% compared to ~ 11% as calculated by garfield++ using the microscopic MC method.
I think even without the exact calculation of the electric field taking into account the actual mesh geometry, is quite good approximation. I do not know if you would like to commend about… I was not aware that AvalancheMC class could do the same with Microscopic NC but much faster… I would like to have a look.
So, I could keep the structure of the Microscopic MC method input file and change only the plate from pad to mesh from MicrocopicMC to AvalanceMC?

Best regards

Georgios

sorry I did not quite understand what you mean with this phrase.

In your last plot, what are the red and blue lines? red = data and blue = simulation with AvalancheMC? Is the agreement with data good enough for you?
greets
Piet

1 Like

Hi Piet,

The micromegas detector we are using is a segmented one. In the last plot is the experimental amplitude spectra obtained with a Fe55 source as we recorded from the mesh and pad superimposed… both spectra should be ideal at perfect construction…
At the best the experimental spectra should give a resolution ~ 14% as our preliminary analysis indicate, compared to ~ 11% we obtain using the microscopic MC class of Garfield++ with this simplified approximation of the mesh used in accordance with your suggestions…
At further stage we need to implement the mesh actual geometry → exact EF evaluation.
I was using NeBEM for that with the standalone application of garfield at CERN but as I reported no access since a couple of months at the executable in lxplus. Unfortunately there is no way to run it at my personal laptop due to the lack of the CERNLIB…
Any solution to keep on using the standalone garfield?

Thanks so much in advance

Best regards

Georgios

Hi,
I hadn’t realised that classic/Fortran Garfield is not available on lxplus at the moment; we should indeed fix that.
You can also use neBEM through Garfield++ though.

1 Like

Dear expert,

Thanks so much. It would be great if this could be fixed since most of my code inputs and results are obtained using the standalone garfield application.

Best regards

Georgios