gemc / source Goto Github PK

Can we add an option to gemc to define the primary particle by kinetic energy instead of momentum? And specify the spread in energy as well. This would be a great help for the momentum-impaired.

Would also like to be able to specify an input spectrum to select the primary from but this would be taken care of if the general particle source is fully implemented so that feature can wait for the GPS.

Can you add an option for Vacuum with optical properties in the HALL_MATERIAL gemc option? It would be just like the Air_Opt option but without the air.

Either add the process if possible, or at least check that the particle is a generated particle and assign procID != 90.

In the code, if NP is less then NBUNCHES, you do not get a LUMI beam at all. Also, if NP is not an integer number times NBUNCHES, all the extra particles go into the last bunch. This is not good for smaller statistics beams. The non-integer extra number of particles should be randomly distributed over the bunches.
There are several solutions:

1. Throw a random number between 0 and 1, and if this number is less than the remnant fraction, add an extra particle in the bunch. I chose this method, since it is closest to the physics. It can cause the number of particles in the event to not equal NP, but on average will give the correct luminosity
2. Before the loop, fill an array with which bunches should get the extra particle, randomly chosen. Then check for each bunch if it needs to get the extra particle. This would guarantee a correct NP in the event, but is more cumbersome to program.
3. Start with method 1, but then adjust the probability as you go from bunch to bunch, so that you always end up with with exactly the correct number of particles in each event. This can be done so the bunches are +/- 1 particle, or it can be done so that each bunch has correct statistical fluctuations in each bunch, but now with the probability changing in such a way that the total is NP over the event.

The request is to implement something like that last option.

Would it be possible to assign magnetic field to a single cad instead of using HALL_FIELD for cad based simulations? This could be assigned as a property for the volume, using the cad.gxml file.

I notice gemc won't quit if it encounter a wrong pid in input
it only give a warning like "Beam Settings >> Particle id 210 not found in G4 table."
Maybe quitting is a safer way?

see http://hepmc.web.cern.ch/hepmc/releases/HepMC2_user_manual.pd

In identifiers.cc, if vector get_identifiers(string var) is passed a string that doesn’t conform to what it expects, it can return garbage. In the particular case I was looking at var was “no”. Then in the loop iden.rule and iden.id are not set by the >> operator and retain an uninitialized value, which is returned in the result.

Arguably this function should print an error message if var is not of the right format, but at the least iden should be initialized.

here is how to repeat this for anyone on ifarm

source /group/clas12/gemc/environment.csh 4a.2.3
gemc /group/clas12/gemc/4a.2.3/clas12NF.gcard -INPUT_GEN_FILE="LUND,/work/clas12/zwzhao/software-tools/local/jpsi_gen_qf_1.txt" -OUTPUT="evio, jpsi_gen_qf_1.evio" -RUNNO=10 -USE_GUI=0 -N=100 -FASTMCMODE=10

it's the same for 1, 2, and 20 mode

We need a detailed info from LUND files (includes info on decaying particles and their parents that are not tracked by gemc) in the final banks, so we can access them later in the analysis stage.

If this is supposed to be the "official" repository for gemc, used by the clas12-offline-software, it should be moved over to the jeffersonlab github so it can be worked on by everyone.

currently gemc searches the default and named field dir and loop through all file "*.dat".
It's not needed this way, instead it can just look for the field files requested.
Also, if there are a few field files have same field name inside, gemc will always use the last one.
This creates unnecessary confusion.

sometimes local field maps are not convenient because it needs a volume.
for example, in EIC, a beamline magnet can have a volume with a local map including empty space in middle where a beam pipe needs to passing through.
Then this beampipe would overlap with the magnet volume.
The only way to avoid this now is to make beampipe in pieces when it go through many magnets and this becomes tedious.
An easy solution would be make those beamline magnet maps all global maps.

If we can have multiple global field map and stitch them together, it would be great for many things.

optMap["ROTATE_FIELDMAP"].help += " Example: -DISPLACE_FIELDMAP="srr-solenoid, 0deg, 0deg, 2*deg"\n";

in file
src/gemc_options.cc

activating optical photons and using the HTCC clas12 detector cause code hangs.
new geant4 patch may fix this

The following bit of code in MPrimaryGeneratorAction.cc is, like, a really bad idea. Introduced on 5/12/2016, so no one checked since then, or using a beam that is off-center is quite rare.
The following loop (line 747):

	do {
		lvx = L_vx - L_dvr + 2*G4UniformRand()*L_dvr;
		lvy = L_vy - L_dvr + 2*G4UniformRand()*L_dvr;
	} while (sqrt(lvx*lvx + lvy*lvy) > L_dvr);

First, it is horrible to spend potentially a massive amount of CPU time on getting an ellipse filled. I.e. if L_dvr is small, but worse, if L_vx or L_vy are larger than L_dvr, this is never satisfied and you enter an infinite loop.That will always be the case if L_dvr is zero, which ought to be a valid number, i.e. no spread.

The correct method for filling a uniform circle is to throw a random radius and a random angle between [0,2*pi]. This is what used to be the case. However, to get a UNIFORM distribution over the circle, you need to make the radius size 1, then take the square root (to get uniform distribution over circle size 1) and then scale up to the desired size. It would be even better if there was an L_dvr_x and L_dvr_y for elliptical beams, which is what we actually usually have in the hall.

hitCollection->insert(noiseHits[h]);

seems not to add hits that integrateRaw process.

GEANT4 seems to be moving toward the GPS over the /particle/gun. The GPS supports adding multiple sources with probabilities so that for any given run, the primary particle is drawn from the distribution of sources. It is possible that way, though a little awkward, to build up a solar system abundance source that looks much like the ions in solar energetic particle events or galactic cosmic rays. Would it be possible to add a feature like this in a future version? Or maybe there already is a way to draw primaries from a distribution of possible initial particles?

feature request, not a high priority

replicas are not working when a volume is placed inside replicas - see attached files.
an improvement of the api would be to add the replica pars in the type description, this will remove the need of an additional volume.

replica.gcard.txt
solid_PVDIS_ec_forwardangle__geometry_Original.txt

make it possible to add Birk constant (https://en.wikipedia.org/wiki/Birks%27_Law) for any external defined material so that its energy deposition is correct no matter optical for scintillation process is turned on or not.

something like

SPREAD_V="even,dvx,dvy,dvz"
SPREAD_V="gaus,dvx,dvy,dvz"

gaus distribution can be obtained easily with "G4RandGauss::shoot(vx/mm,dvx/mm)"

At some point would be nice to add electric fields along with magnetic. Electrostatic fine to start with, changing fields farther off?

It would be nice to have the option to shift and rotate field. This allow simple field map in some cases

If use pi+ as the beam, K+ and p as background, generated vertex distribution of K+ and p are not the same as input one. In gcard file, I have following commands:

But the generated vertex for K+ and p are much wider than input. Please find an example plot in attachment.
vertex_rich_pikpplus_flat.pdf

It would be nice to add a feature to calculate RadiationLength and NuclearInterLength along certain path.
Maybe shoot geantino and use G4Material::ComputeRadiationLength and G4Material::ComputeNuclearInterLength() to sum the materiel along its path.
I wonder if checking overlap is done in a similar way. If so, it can be expanded easily.

Please add support for 3D magnetic field map, it can be in XYZ or ZRPhi coordinate. maybe RThetaPhi coordinate is not really needed.

Hello,

For the Micromegas (BMT, FMT, FTM) as well as the RTPC we would need to be able to fill multiple ADC as a function of the time for each event. At the END of one event, we would integrate the signal left by all the hits that left a signal on one strip. In other words:

Event 42:
Hit 1: would leave on pad N ADC_1(t)
Hit 2: pad N, ADC_2(t)
Hit 3: pad N, ADC_3(t)
...

And at the end of the event, pad N would see sum(pad N, ADC_X(t))

Where t and N are of course adjustable parameters. I suppose I can take care of the sum part in hitprocess but I need to be able to output ADC(t) for each hit with the shape function of my choice.

This step is very important for us to evaluate precisely occupancy, resolutions, fitting procedures...

Thank you!

in file "source/materials/material_factory.cc"

// helium 3 material gas
// Density at 21.1°C (70°F) : 0.1650 kg/m3
mats["helium3Gas"] = new G4Material("helium3Gas",  0.1650*mg/cm3, 1, kStateGas, 294.25*kelvin);
mats["helium3Gas"]->AddElement(helium3, 1);

The density seems a bit large
0.932_44.6(amg=mol/m3)_3.016(g/mol)=1.254e-3g/cm3

refer to
https://en.wikipedia.org/wiki/Amagat

Would like to be able to select only certain values in the true information bank to save. Many times one only needs one or two variables. Could apply to the digitized banks also?

When gemc is built with newer versions of gcc, with c++11 support, there is a runtime error.

Entering initializeCNDConstants
Entering initializeCTOF
  > Data loaded in translation table dcTT
 > bmt pre-initizialization.
terminate called after throwing an instance of 'std::bad_alloc'
  what():  std::bad_alloc

The error first appears when initializing the global object bmtConstants in source/hitprocess/clas12/micromegas/BMT_hitprocess.cc (line 332).
bmtConstants BMT_HitProcess::bmtc = initializeBMTConstants(-1);

It seems to be a memory allocation issue and may also be an issue with older versions of gcc but before stdc++11 the error was not reported.

I would like to request the option to not specify the number of simulated event in a LUND file. The LUND file could be loaded as usual and events would be simulated until the EOF is reached. This would be very helpful in the event you are given a generator LUND file and the total number of events is not specified.
Thanks!

after ./go_all

show:

!! Attention: GEMC installation check /home/inat/APP/jlab_software/devel/Linux__Ubuntu16.04-x86_64-gcc5.4.0/gemc/devel/gemc not found

node1:~/APP/jlab_software/devel/Linux__Ubuntu16.04-x86_64-gcc5.4.0/gemc/devel> ls

api futureDevelopments materials README.md build_log gemc.cc mirrors SConstruct changeLog generator output sensitivity detector gui parameters src fields hitprocess physics utilities
`

What is my wrong？

We can use 1/8 of in XYZ 3D to represent a dipole field map with reflection symmetry. namely, (0,xmax)(0,ymax),(0,zmax) can represent (-xmax, xmax)(-ymax,ymax),(-zmax,zmax)

I recommend using this for every sensitive logical volume. This would allow for more flexible hit processing and output.

From http://geant4.cern.ch/support/ReleaseNotes4.10.2.html

Digitization & Hits
Added new G4MultiSensitiveDetector functionality, allowing to assign multiple sensitive-detectors to a single logical-volume. Sensitive-detectors are added to this proxy class and an instance of the proxy is assigned to the logical volume. Calls to the sensitive-detectors methods are forwarded to all user-defined sensitive-detectors that are added.

I could be wrong, but given that
SNR=NPE/sigma, then sigma = NPE/SNR. With the RCA handbook giving SNR=sqrt(NPE)/pmtFactor, shouldn't the code linked above be calculating sigma as sigma=sqrt(NPE)*pmtFactor?

In the source it mentions

Note: GSIM sigma was incorrect (used 1/sigma for sigma).

But I am not sure quite what this refers to.

a feature request: changing material at run time.
This could avoid a lot of duplicate geometry file and perl scripts.
I think it can be done in two ways

1. by giving a volume name, change its material to a new material
2. by giving a material name, change all volume with that material to a material.
   both can be useful at different situation.

by default, particle with status 0 in lund won't be used in simulation, nor included in output "generated" bank. It would be nice to have an option to just include it output "generated". This is useful for example, if we want to pass info of target nucleon inside nucleus before scattering into output

Have a special bank that can record the calibration constants onto the data stream.

Several applications and studies require knowledge of particle decay info (daughter vectors of primary particle decays):

1. Displaced vertex studies.
2. Charged pion decay in CLAS12 volume.
3. Neutral pion decay for two photon reconstruction in calorimeters.

Perhaps a special bank that keeps track of primary particle decay kinematics could be added.

Cole Smith

I would like to be able to generate a specific ion (Fe-56 or C-12 or whatever) as a primary particle in gemc. There is a GenericIon in the particle list but need to select what ion.
I have tried to set the ion in the G4dialog tab using:
/gun/particle ion
/gun/ion 26 56 2
/gun/position 0 0 2 cm
/gun/direction 0 0 -1
/gun/energy 1 GeV
/run/beamOn 1

I may not have the verbosity set right but have seen no indication that any message was reported with any of these commands. However, if I forget to set '/gun/particle ion' first I do get an error message when I try '/gun/ion 26 56 2'. I've also tried it without the '2' charge state. Same result if I use the gemc 'Run' button after setting these instead of '/run/beamOn 1'.

I tried this again but setting the main particle to GenericIon in the gemc drop down menu. I do see messages showing the particle changed to GenericIon in the output but nothing further when I try to set the ion using then G4dialog commands.

I'm using a simple time-of-flight as an example and for a fixed energy, the times between detector hits do not appear to change depending on what mass ion I'm trying to set.

I'm having trouble setting up a magnetic field using the dipole or cylindrical symmetry options. Separately, I am also not able to load a magnetic field map. I think it is very likely that I'm just doing something wrong but am running out of guesses. Can you help? I also don't fully understand the parameters and possible values available to use in the field definition file. Sometimes I see bx, by, bz, units, and other times Bunit, etc. This would be a good place to add some documentation if the mood strikes... :)

1. simple fields
   I have a small example that creates a field along the z-axis using 4 different definitions in separate volumes. I have a uniform field, a multipole (n=2) field rotated to point along z, a multipole (n=4) field also rotated, and a solenoid field defined using the cylinder-z symmetry. With the first 3 definitions I get a field created and can see the arrows with the camera button to show magnetic field. The solenoid field says it is created in the gemc output but no field is shown in the camera option and electrons don't bend in the volume. Any idea what I'm doing wrong here? Example is in dipole_magnet.zip

2. field maps
   In this example I am trying to define a field using a field map in the field definition file. I wrote a small python executable to generate the field file which should be aligned with the z-axis and uniform across a 20cm x 6cm x-y plane except that it reverses direction for x<0. The example finds the field definition and reads the map values correctly. It builds the field, attaches it to the volume. Then when a particle tracks through the volume gemc reports that is tracking in the field, but with large uniform steps and reports the field as 0. Using the camera option to look at the field gives an error message that there is no field in the scene. This may be related to the symmetry, not sure. The example is in ex2_field_map.zip

thanks!

dipole_magnet.zip
ex2_field_map.zip

I could really use the "gflux" sensitive detector type we've discussed that would count numbers of particles entering the volume without need to track in the volume. This would be used for something like a PMT or SiPM where there may be many thousands of optical photons arriving at different points and directions but all I care about is how many were detected. It should be general for all particles, though, not just for optical photons.

A variation on this would be to save the time structure, for example, numbers of photons vs. time in the event. This might be done by treating Nphoton vs time like a voltage waveform and then reading it out through the same way as the FADC. Just one possibility. The benefit would be for something like a scintillator with a Cerenkov component - one could sample the fast and slow components without having to save hundreds or thousands of hits, in some cases. Or is there a better way to do that?

For many purposes it can be useful to be able to re-run the simulation for a small subset of events. For example, you might have a large number of simulated events out of which there are ten which you'd like to examine visualization of, or which you'd like to re-run with mother saving turned on.

For SoLID at one time (several years ago) I had a kludged version of this working. There were command line options to save the simulation state to a file, and to re-run the simulation for only the events whose state had been saved. The way saving worked was that before each event it would internally store the state of the RNG and (if any) the current input file data. During the event, the hit processor could call a routine called "SaveMe". After the event, if "SaveMe" had been called, the RNG and input data would be written to their own output files. To replay, if the command line option specified it, the RNG would be set on each event from the next saved state and the primary generator would be initialized from the next saved input data.

This worked for my purposes but was rather a kludge, since if you wanted to change the criteria for saving events you had to change the hit processor code and recompile.

Two ways to avoid that problem would be: (1) Provide a way to specify a selection criterion as an input, rather than hard coding it, or (2) Save the RNG for all events, say in the header bank or in its own new bank; then to replay, specify numerically which events should be re-read from the input file and have the RNG set from the stored information in the output file.

These have obvious drawbacks: (1) would require a parser similar to what ROOT uses for tree selection, a major piece of work. (2) requires adding, I believe, 625 unsigneds (the size of the RNG state) to the output for every event (if requested), a high storage price to pay. It may be the best way to do this, though. One presumably would use this option infrequently and not on enormous data sets, and 625 unsigneds per event would typically not be a large fractional increase in the size.

However, it would then be necessary to provide the capability for GEMC to read its own output files, potentially another can of worms.

Should this condition every be true? The mirror handling is done in sensitivity not hitType, so this should never happen.

It seems minStep,integration,shift are not supported in map file file any more. They are simply ignored if defined in a map file
All of them have to be set in option.
This break backward compatibility at least from v2.3 when they are not ignored in map files.
And this can make gcard file unnecessary long

I would suggest to keep support of minStep,integration,shift, and add rotation in map files.
Then make options to override whatever defined in map.
This would restore backward compatibility and provide most flexibility.

I can try to implement them. Let me know.

Just a reminder for version 3.x ---
really want to eliminate the requirement in 2.x versions that gemc has to resolve the jlab MySQL host in order to run even if it isn't using MySQL or anything at JLab.

Just ran into that problem again and thought I'd put it on the list so we don't lose track of it.

thanks!

modify source/physics/PhysicsList.cc to allow hadron physics list only without EM physics list to run

This can used to test geant4 physics

Currently -N=100 will run 100 events even when input file only has 50 events
another feature would be nice
When using input file, gemc don't need -N option and it can will run all events in the input by default
Add N option only if you want run number of events smaller than input file.
It will never run number of events larger than input file

Currently the internal particle gun doesn't shoot evenly in solid angle, it would be nice to do so

this will do the trick
change this line in src/MPrimaryGeneratorAction.cc
Theta = stheta/rad + (2.0_G4UniformRand()-1.0)_dtheta/rad;
to
Theta = acos(G4UniformRand()*(cos(theta/rad-dtheta/rad)-cos(theta/rad+dtheta/rad))+cos(theta/rad+dtheta/rad))/rad;

will need to check if this works with option ALIGN_ZAXIS