NPONT Error for High Reduced Electric Field Gas File with Magboltz

[NSos]

Hi! I am currently working on calculating gas files with electric fields in a thin isobutane gas with thin wires at a high electric field for electron multiplication purposes. I get up to around 20-30 kV/cm, but then get an error

PROGRAM STOPPED NPONT= 8001 ITER= 379073975

I tried modifying magboltz.f to have larger arrays than 8000 and making the warning come up for larger values and recompiling Garfield (a desperate attempt at addressing this, I must admit), but the warning just comes up for whatever higher number I set. Any help in troubleshooting this would be appreciated!

The full output looks like this:

Processing 99 wires along x and y axis each.

   Starting wire position -4.9
MediumMagboltz::SetComposition: iC4H10



Sorting out the gas table


Writing gas file for E = 30000 V

MediumMagboltz::GenerateGasTable: Found 15 excitations and 3 ionisations.
    iC4H10 EXC.  TRIPLET                    ELOSS=   7.2    , energy = 7.19993 eV.
    iC4H10 EXC.  TRIPLET                    ELOSS=   8.6    , energy = 8.59992 eV.
    iC4H10 EXC. SINGLET         F=.0013     ELOSS=   7.5    , energy = 7.49993 eV.
    iC4H10 EXC. SINGLET         F=.0150     ELOSS=   8.0    , energy = 7.99992 eV.
    iC4H10 EXC. SINGLET         F=.1140     ELOSS=   8.5    , energy = 8.49992 eV.
    iC4H10 EXC. SINGLET         F=.1570     ELOSS=   9.0    , energy = 8.99992 eV.
    iC4H10 EXC. SINGLET         F=.1710     ELOSS=   9.5    , energy = 9.49991 eV.
    iC4H10 EXC. SINGLET         F=.1880     ELOSS=  10.0    , energy = 9.99991 eV.
    iC4H10 EXC. SINGLET         F=.2050     ELOSS=  10.5    , energy = 10.4999 eV.
    iC4H10 EXC. SINGLET         F=.1930     ELOSS=  11.0    , energy = 10.9999 eV.
    iC4H10 EXC. SINGLET         F=.1620     ELOSS=  11.5    , energy = 11.4999 eV.
    iC4H10 EXC. SINGLET         F=.1030     ELOSS=  12.0    , energy = 11.9999 eV.
    iC4H10 EXC. SINGLET         F=.0670     ELOSS=  12.5    , energy = 12.4999 eV.
    iC4H10 EXC. SINGLET         F=.0640     ELOSS=  13.0    , energy = 12.9999 eV.
    iC4H10 EXC. SINGLET         F=.0280     ELOSS=  13.5    , energy = 13.4999 eV.
    iC4H10 IONISATION                       ELOSS=  10.67   , energy = 10.6699 eV.
    iC4H10 IONISATION-EXCITATION (BREAKUP)  ELOSS=  17.0    , energy = 16.9998 eV.
    iC4H10 IONISATION  K-SHELL              ELOSS= 285.0    , energy = 284.997 eV.
MediumMagboltz::GenerateGasTable: E = 30000 V/cm, B = 0 T, angle: 1.5708 rad
 RM48 INITIALIZED:  54217137           0         0


          PROGRAM MAGBOLTZ 2 VERSION 11.19

          MONTE CARLO SOLUTION FOR MIXTURE OF  1 GASES.
     ------------------------------------------------------

       GASES  USED                 PERCENTAGE USED 

      ISO-C4H10 2014 ANISOTRPIC      100.0000


  GAS TEMPERATURE =  20.0 DEGREES CENTIGRADE.
  GAS PRESSURE =     5.3 TORR.

  INTEGRATION FROM 0.0 TO ******** EV.  IN 4000 STEPS. 

 PENNING EFFECTS NOT INCLUDED

 ANISOTROPIC SCATTERING TYPE 2 (OKHRIMOVSKYY) USED IF AVAILABLE

 SHORT DECORRELATION LENGTH = 400000 COLLISIONS.


 THERMAL MOTION OF GAS INCLUDED

  ELECTRIC FIELD =  30000.0000 VOLTS/CM.
  MAGNETIC FIELD =     0.0000 KILOGAUSS.
  ANGLE BETWEEN ELECTRIC AND MAGNETIC FIELD =    90.000 DEGREES.
  CYCLOTRON FREQ. =   0.000D+00 RADIANS/PICOSECOND

  INITIAL ELECTRON ENERGY =5242.880 EV.

  TOTAL NUMBER OF REAL COLLISIONS =   40000000

  NULL COLLISION FREQUENCY FOR EACH GAS COMPONENT IN UNITS OF  (*10**12/SEC)
         0.242D+00         0.000D+00         0.000D+00
         0.000D+00         0.000D+00         0.000D+00


    VEL      POS        TIME      ENERGY   COUNT   DIFXX     DIFYY     DIFZZ

 10854.89  0.216D+03  0.199D+08 953.7614    10   ********  ********       0.0
 10878.53  0.434D+03  0.399D+08 959.2296    20   ********  ********       0.0
 10842.45  0.649D+03  0.599D+08 948.3591    30   ********  ********  ********
 10840.89  0.865D+03  0.798D+08 947.1731    40   ********  ********  ********
 10832.42  0.108D+04  0.997D+08 945.6626    50   ********  ********  ********
 10831.58  0.130D+04  0.120D+09 944.2754    60   ********  ********  ********
 10835.47  0.151D+04  0.140D+09 945.1777    70   ********  ********  ********
 10835.25  0.173D+04  0.160D+09 944.7176    80   ********  ********  ********
 10838.04  0.195D+04  0.180D+09 946.1796    90   ********  ********  ********
 10833.12  0.216D+04  0.199D+09 944.9686   100   ********  ********  ********
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------

  CALCULATED MAX. COLLISION TIME = 278.02 PICOSECONDS.

  NUMBER OF NULL COLLISIONS =    8316124
  NUMBER OF REAL COLLISIONS =   40000000

  Z DRIFT VELOCITY = 0.1083E+05 MICRONS/NANOSECOND  +-    0.11% 
  Y DRIFT VELOCITY = 0.0000E+00 MICRONS/NANOSECOND  +-    0.00%
  X DRIFT VELOCITY = 0.0000E+00 MICRONS/NANOSECOND  +-    0.00%


           DIFFUSION IN CM**2/SEC.


  TRANSVERSE DIFFUSION   = 0.7602D+07 +-    9.31%
          =210.53299 EV. +-   9.308%
          = 1184.717 MICRONS/CENTIMETER**0.5  +-    4.65%


  LONGITUDINAL DIFFUSION = 0.2023D+09 +-    15.0%
          =5601.3565 EV. +-   14.96%
          = 6110.841 MICRONS/CENTIMETER**0.5  +-    7.48%



  IONISATION RATE /CM.= 0.8309D+02 +/-  0.02 PERCENT.
  ATTACHMENT RATE /CM.= 0.0000D+00 +/-  0.00 PERCENT.



  MEAN ELECTRON ENERGY = 944.9686 EV. ERROR =  +-    0.39%


  SOLUTION FOR STEADY STATE TOWNSEND PARAMETERS
  -------------------------------------------------

 SPACE STEP BETWEEN SAMPLING PLANES = 0.15556D+03 MICRONS.

   PROGRAM STOPPED NPONT= 8001 ITER= 379073975

And the relevant code is here:

	const float P_mbar = 7; //7 mbar isobutane in the volume	
	const float P = 0.7500617 * P_mbar; //Pressure needed in in torr
	const float V_wire = 650.;
	const double radius = 0.003; //Wire radius
	const double halflength = 5.; //Wire half-length
	const double chamber_depth = 1.;
	const float pitch = 0.1; //Distance between wires in cm
	//const float mylar_thickness = 0.00009;
	const float mylar_thickness = 0.;
	const float plot_z_height = 0.6;
	const float wires1_height = 0.3;
	const float wires2_height = 0.9;
	const float cathode_height = 0.6;
	
	//Geometry
	GeometrySimple geo;

	std::vector< std::vector<SolidWire> > Wires; //Vector to list wires

	Wires.resize(2); //Column 0 is for x wires, column 1 is for y wires

	float start_wire_pos = -1. * halflength; //Position in cm of first wire in cm in along +/-5 cm axis
	float end_wire_pos = halflength; //Same as above but last wire

	int N_wires = floor((end_wire_pos - start_wire_pos)) / pitch; //Number of wires fitting into that space

	MediumConductor metal;
  
  //Iteratively calculate x and y wire plane positions and add them to vector
	for(int i = 0; i <= N_wires; i++){
  	
		float wire_pos = start_wire_pos + (float)i * pitch;
  	
		if(i == 1){std::cout << "   Starting wire position " << wire_pos << std::endl;}
  	
		SolidWire wire1(0., wire_pos, wires1_height, radius, halflength, 1, 0, 0);
		SolidWire wire2(wire_pos, 0., wires2_height, radius, halflength, 0, 1, 0);
		
		wire1.SetBoundaryPotential(V_wire);
		wire2.SetBoundaryPotential(V_wire);
		
		Wires[0].push_back(wire1);
		Wires[1].push_back(wire2);
	}
  
	//Add wires to geometry
	for(int i = 0; i < (int)Wires[0].size(); i++){
		
		geo.AddSolid(&Wires[0][i], &metal);
		geo.AddSolid(&Wires[1][i], &metal);
	}

	//MediumConductor mylar;
	MediumMylar mylar;

	//Create a solid cathode plane made of Mylar between x and y wire planes
	SolidBox cathode(0., 0., cathode_height, halflength, halflength, mylar_thickness / 2.);

	//Cathode is grounded
	cathode.SetBoundaryPotential(0.);

	geo.AddSolid(&cathode, &mylar);

	//Set up gas
	MediumMagboltz gas("ic4h10", 100.);

	gas.SetPressure(P);
	gas.SetTemperature(293.15); //Kelvin
	gas.EnableThermalMotion(true);

	std::vector<double> Es = {30000., 40000., 50000., 100000., 200000., 300000.};
	
	const int ncoll = 1;
	
	std::vector<double> bfields = {0.};
	std::vector<double> angles = {1.570796};
		
	for(int i = 0; i < (int)Es.size(); i++){
		
		std::vector<double> efields = {Es[i]};
			
		std::cout << "Writing gas file for E = " << Es[i] << " V" << std::endl << std::endl;
			
		gas.SetFieldGrid(efields, bfields, angles);
		gas.GenerateGasTable(ncoll);
			
		std::string gas_fname = "Gas_Tables/isobutane_";
		std::stringstream nameformat;
		nameformat << Es[i];
		gas_fname += nameformat.str();
		nameformat.str("");
		gas_fname += ".gas";
			
		gas.WriteGasFile(gas_fname);
			
		std::cout << std::endl << std::endl << "Written" << std::endl << std::endl << std::endl;
	}
	
	geo.SetMedium(&gas);

[hschindl]

Hi,
30 kV/cm at a pressure of 7 mbar is a very high field. Are you sure that you have such a high field in your detector?

[NSos]

I plot reduced field in the z-axis direction elsewhere in my code, attached below, and it suggests much higher fields than that in the vicinity of the wires, as the simulation I am running is an MWPC. Furthermore, elsehwere in my code, I have AvalancheMC calculation, as I am trying to calculate electron current multiplication. For that reason, I was hoping to have some fairly reasonable estimates of gas properties even in high fields in the vicinity fo the wires, so a few hundred of kV/cm.

At the moment, for electrons near the wires, I get errors from AvalancheMC, where it cannot calculate velocity, and I think that’s because my highest gas table reduced field is 20 kV at the moment. If you think I am misusing Garfield for such high fields, I can probably make do with what I already have, but I was hoping to explore this in some more detail.

[Piet]

Hi @hschindl

I was trying to see wether with ViewMedium one could read data 100V/cm-20kV/cm from the gasfile and then make the plot with an exponential extrapolation outside the 20kV/cm range, but that was not possible. Asking from the Medium to have swarm parameters for > 20kV/cm was not possible. I believe that in this case that would be useful as it would be a first reasonable estimate for the swarmparameters in the field close by the wire. What is your thought on this?
greets
Piet

[hschindl]

Hi,
can you try running Magboltz without the “gas thermal motion” option (i. e. commenting the line gas.EnableThermalMotion(true);)?

Btw, is it really the reduced field that’s shown in your plot?

[hschindl]

Hi Piet,
hm, that’s strange. If you ask for a field outside the range of the table, Medium should always try to do an extrapolation. If there are at least two fields in the table, that should normally work. Are there any messages indicating why the extrapolation failed?

[NSos]

I make that plot using command

ViewField fieldView_E(&nebem);
fieldView_E.SetNumberOfSamples1d(1000);
fieldView_E.PlotProfile(-1. * pitch, 0., 0., pitch, 0., chamber_depth, "E", 0);

Here, pitch is 0.1 cm wire separation, and chamber_depth is 1 cm depth of the detector volume in the z-direction. The plot I sent is simply the output of the command here, I did not adjust it in any way. I called it reduced field, as its units are listed as V/cm. I appreciate “reduced field” can also mean Volts per unit pressure as well, but I do not know how to make Garfield calculate one of them, and I don’t particularly need it. My understadning is, however, that gas tables calculated Magboltz are using this V/cm parameter specifically, which is why I wanted to calculate it to greater “reduced” fields and encountered the error above.

EDIT: Running the gas file calculation now without thermal motion. Will update once it’s done.

[NSos]

Just as an update, the code without thermal motion has been running for 6 hours, but so far has not calculated a gas file for 30 kV. However, it also would have produced the error I mentioned a long time ago, so something is definitely working better. For reference, I run single-core on a 13th Gen Intel Core i7-13700H CPU. The process is definitely still going; if I look at top output, it’s listed as using 100% of 1 core.

[hschindl]

Hi,
if it’s still running that’s in some sense promising, but such a long run time probably also means that you have a huge gain/Townsend coefficient. I’m still puzzled by the high reduced electric field in your detector. Is this really the voltage configuration at which you operate the detector?

[Piet]

Hi all

I called it reduced field, as its units are listed as V/cm. I appreciate “reduced field” can also mean Volts per unit pressure as well

@NSos please do not create confusion about definitions and/or units. There is no universe where we will call a normal electric field (in V/cm) a “reduced” electric field as this term is reserved for something else. Please have a look here.

Potential (V) => units: Volt (V)
Electric Field (E) => units: V/cm
Reduced Electric Field (E/N) => units: Td

@hschindl , I think that @NSos showed a plot of the electric field, which is reasonable for a field close by the wire.

@hschindl: For what concerned my question, I do not recall exactly what I was doing. Maybe I was asking higher values to the fieldgrid, which of course is not possible. When using the ElectronTownsend function of class Medium, then it works fine to obtain an estimated Townsed value outside the range of the fieldgrid provided by the gasfile. However I was not able to use the ViewMedium class plotting tools, changing the range by view.SetRangeE() or view.SetRangeY(), I had to do it manually. I tried to use SetRange both before and after the Plot function, but it did not make any difference. I attach a small piece of code that reproduces my problem. Maybe I do not know how to use it correctly?

@hschindl I attach my extrapolation for the Townsend coefficient, but I have my doubts this is ok. I used linear extrapolation, but I would think we should add some “logaritmic” extrapolation option (right now available: constant, linear, exponential). What’s your thought about this?

update: having linear scaling is likely good enough, however as you can see from this plot this scaling is not linear.

reproduce.C (631 Bytes)

[hschindl]

Hi @piet,
to check the scaling can you plot α/p as function of E/p?

The “exponential” extrapolation is in fact an extrapolation on a log-scale:

Source/Medium.cc · master · garfield / garfieldpp · GitLab

Cheers,
H.

[Piet]

Thanks Heinrich for the suggestion.
Here below you can see α/N vs E and α/N vs E/N.
Hope plotting vs /N instead of /p is fine for you?

I read the code in Source/Medium.cc (lines 1487-1514) [1], to me it seems there is nothing obviously wrong, however when you look at the log-log plot, you can see that the linear extrapolation is curved (in a way similar as the exponential extrapolation). A linear function plotted in a log-log plot should be a straight line, which is not the case here.

Any ideas? Thanks!
Piet
[1] Source/Medium.cc · master · garfield / garfieldpp · GitLab
Piet

[hschindl]

I see what you mean… I think that’s because for the Townsend coefficient α (and also the attachment coefficient η), the gas table contains log(α), so it’s a log-log interpolation of log(α)…

[Piet]

Thanks @hschindl, I needed some time to digest. So it is a linear interpolation using log(α) values saved in the gas table. I made a fix for it in the Medium class, and this is the result, which is much better, but maybe not yet perfect. If you are ok I make a merge request with this feature, you can have a look.

@NSos, for what concerns your question, I think you can calculate the gasfile up to 20kV/cm (at 7mbar) and use the extrapolation of Townsend coefficient to higher fields.
In order to know whether the simulation provides reasonable values, it would be good (not to say necessary) to verify the actual gain in measurements. Do you have data to compare to?

greets
Piet

[NSos]

Hi @Piet, thank you very much for all the checks. We do have data, that’s being collected by students now, so in a few weeks, I can probably do a comparison.

One of the reasons I wanted to have a look at these high fields is because I get errors getting velocity from AvalancheMC near the wires. Another issue there is that when the electron cloud grows the programme slows down a lot (unsurprisingly). So at the moment the programme runs too slow for me to find out what the multiplication is like for electron produced further from the wires, as they are still drifting without multiplying by the point the closer electrons start blowing up the processing time. I was hoping better gas tables may accelerate this.

If I get around this (or just accept waiting for days/weeks for processing), I was planning to calculate total drift time of all the electrons, and then add signal calculation to my sensor, so that I can compare the pulse shape.

Do the changes in processing time I describe make sense to you and match what you would expect? If you wish to check it, I can provide the source code.

Another question, if you are merging some changes into your github, would you suggest redoing my previous gas tables after cloning them to my machine? Or lower-voltage tables would be unaffected?

[Piet]

Dear @NSos

• If AvalancheMC is running to slow, in general, I wouls suggest you try driftLineRKF.
• Slow simulation can also be because electrons are trapped in a field with nearly 0 Efield,
  so it would be good to verify your electric field map.
• It is always good to have a look at the source code to see if there is no suboptimal approach used
• Previous gastables are OK, as the modifications are only targeting the extrapolation methods used on the data in the gasfile.

kind regards
Piet

[NSos]

OK, @Piet , that’s great, thanks for the explanation! I will try driftLineRKF, hopefully this will solve the matter. And thank you and @hschindl for this investigation!