JeffersonLab / remoll

Simulations for the MOLLER Experiment at Jefferson Lab, http://moller.jlab.org
http://jeffersonlab.github.io/remoll/
11 stars 55 forks source link

Aluminium generator beam energy factor #443

Closed cipriangal closed 3 years ago

cipriangal commented 3 years ago

Environment: (where does this bug occur, have you tried other environments)

Steps to reproduce: (give a step by step account of how to trigger the bug)

  1. the generator for both the quasiElastic (https://github.com/JeffersonLab/remoll/blob/master/src/remollGenAl.cc#L201) and inelastic (https://github.com/JeffersonLab/remoll/blob/master/src/remollGenAl.cc#L152) generators have a beam-Energy factor in the x-section
  2. confirmed this was in the original Qweak generators, but not clear why it is there. The original ref (https://arxiv.org/pdf/1203.2262v2.pdf -- eq3) doesn't have this
leafybillow commented 3 years ago

I think the eq3 in the original ref is energy differential x-sect, of which the dimension is [ubarn/(s.r * energy) ] In the generator, the x-sect should be energy integrated and in [ubarn / s.r.], which is relevant for rate weighting.

hanjie1 commented 3 years ago

Timing a beam-energy factor isn't integrating over E'. The cross section we need as weight should be differential cross section, otherwise why we need to generate random scattered energy.

leafybillow commented 3 years ago

Timing a beam-energy factor isn't integrating over E'. The cross section we need as weight should be differential cross section, otherwise why we need to generate random scattered energy.

  1. In remollGenAl.cc, the beginning of SamplePhysics gets 'beam energy' from a remollVertex. This 'beam energy' is not always the same as 11 GeV and is subtracted by radiation loss in material with a probability encoded in remollBeamTarget.cc :: SampleVertex().
  2. In both quasielastic and inelastic generator, nu, W and Q-squares are calculated based on both eOut and aforementioned 'beam energy'. The reason for a flat sampling of outgoing electron energy is we want an uniform(almost) distributed kinematics , which are needed for form factor functions. This is similar to a scan of x-sect over different kinematics.
  3. With a given set of kinematics fixed, integrating the 'eq 3' over all possible E' from electron energy at rest to primary beam energy (not necessarily 11 GeV) yields the ' (beamE/GeV - electron_mass_c2/GeV) ' term you mentioned. And this gives differential x-sect [barn/s.r], which is later multiplied by phaseFactor[s.r] and converted to an effective x-sect in [barn].
  4. Which variable should be used for rate weighting? It depends on how we calculated rate. In remollPrimaryGeneratorAction.cc :: GeneratePrimaries(G4Event*) needs this x-sect in [barn], beam luminosity and target density to calculates the rate, otherwise the dimension of rate is wrong. Chapter 4 in Peskin's QFT( p.100 ) gives the relation between rate, luminosity and x-section.
hanjie1 commented 3 years ago

The cross section represents when you have a particle incident with beamE, scattered out with eOut at angle theta, what the probability of this process happens. DIS cross section is per Omega per E'. We detect particle at specific E' and theta, so there shouldn't be an E' integration. On the other hand, if you want to do E' integration, it's not integrated over to beamE, it should be integrated over to the elastic E'.

hanjie1 commented 3 years ago

Probably we shouldn't rely on Qweak inelastic generator. If you take a look at "Alinel_crsec_gen", which claims to be inelastic, but the mott scattering cross section they use it's for elastic: G4double sigma_mott = pow(Z_Al0.719982/E_inCTH/(STHSTH),2)/(1+2E_in/M_AlSTHSTH)*10000;

jmammei commented 3 years ago

But the value that is output is the "effectiveXsection" which uses the sigma_mott but multiplies by form factors... it looks about right at quick glance...

In the other emails did "timing" mean "times-ing" or multiplying?  I was a little confused by that...

Juliette

On 2021-03-04 6:49 p.m., Hanjie Liu wrote:

Probably we shouldn't rely on Qweak inelastic generator. If you take a look at "Alinel_crsec_gen", which claims to be inelastic, but the mott scattering cross section they use it's for elastic: G4double sigma_mott = pow(Z_Al/0.719982/E_in/CTH/(STH/STH),2)/(1+2/E_in/M_Al/STH/STH)*10000;

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_JeffersonLab_remoll_issues_443-23issuecomment-2D791061032&d=DwMCaQ&c=CJqEzB1piLOyyvZjb8YUQw&r=w6PhlF4C-oeM-6LmjCQ3mg&m=pDLhBWbRmZwXFT816v8qVOuLLsbYQsoITD2dPc6z0Zc&s=kBUq8xMyBZCfmjIAvpGwp0xvH3Pk5gA92-s2FxqFI64&e=, or unsubscribe https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AB4CGEDMU5PKQOYYB3VPPSLTCATA3ANCNFSM4YTWQS4A&d=DwMCaQ&c=CJqEzB1piLOyyvZjb8YUQw&r=w6PhlF4C-oeM-6LmjCQ3mg&m=pDLhBWbRmZwXFT816v8qVOuLLsbYQsoITD2dPc6z0Zc&s=5FaMmnJpXL5mf8DaF0xM0j1gummEfvHnVphDCzZKOIg&e=.

-- Dr. Juliette Mammei

Associate Professor Dep't of Physics and Astronomy University of Manitoba Winnipeg, MB R3T 2N2 Canada

hanjie1 commented 3 years ago

Our issue is in the inelastic generator. As you said there is a form factor in "Alinel_crsec_gen", the cross section couldn't be for inelastic. For the inelastic cross section, we use structure functions.

Yep..I mean "multiplying"...I just assumed you can use times as "timing"....

The reason that there is a "phaseSpaceFactor" as Tao mentioned in effectiveXsection is to remove the non-uniform solid angle when sampling theta and phi. Here is a talk from Tyler Kutz about g4hrs, which has a similar idea: https://prex.jlab.org/DocDB/0001/000124/001/g4hrs_calculations.pdf

jmammei commented 3 years ago

The angle phase space factor is correct (calculated in line 71, applied in line 102) but there is no energy phase space factor there.  You do need that to get the right weight... the inelastic energy is generated from a flat distribution (calculated in line 125), so there should be a Delta E factor somewhere.  Maybe beamE isn't exactly the right factor, but there needs to be one.

If the distribution is flat, then the phase factor is a simple Delta of the max and min values.  If you draw from a non-uniform distribution, then the phase space factor is more complicated. But you still need one even in the flat case.

The inelastic formula calculates W1 and W2 from the F1 and F2 structure functions (not form factors).  Are both the inelastic and quasielastic terms included in line 201?  Also in line 152? Are they calculating the same thing?  Eq 3 from the source you sent earlier?

In any case I am glad someone is looking at this closely...

Juliette

On 2021-03-04 7:45 p.m., Hanjie Liu wrote:

Our issue is in the inelastic generator. As you said there is a form factor in "Alinel_crsec_gen", the cross section couldn't be for inelastic. For the inelastic cross section, we use structure functions.

Yep..I mean "multiplying"...I just assumed you can use times as "timing"....

The reason that there is a "phaseSpaceFactor" as Tao mentioned in effectiveXsection is to remove the non-uniform solid angle when sampling theta and phi. Here is a talk from Tyler Kutz about g4hrs, which has a similar idea: https://prex.jlab.org/DocDB/0001/000124/001/g4hrs_calculations.pdf https://prex.jlab.org/DocDB/0001/000124/001/g4hrs_calculations.pdf

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_JeffersonLab_remoll_issues_443-23issuecomment-2D791081418&d=DwMCaQ&c=CJqEzB1piLOyyvZjb8YUQw&r=w6PhlF4C-oeM-6LmjCQ3mg&m=CgBqx82iG1iHV7o8fEL5ezRMkxxacPKsNcXZ3BlwMwE&s=dzdO-th9KnBg6peh0_e3p8PO9v6PLoeuD2LTizc9fHE&e=, or unsubscribe https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AB4CGEGAD26YFCD3T2B5Q4DTCAZSTANCNFSM4YTWQS4A&d=DwMCaQ&c=CJqEzB1piLOyyvZjb8YUQw&r=w6PhlF4C-oeM-6LmjCQ3mg&m=CgBqx82iG1iHV7o8fEL5ezRMkxxacPKsNcXZ3BlwMwE&s=NQcddGN3oFUHsbt0MdFV3V9WlGnHFzSqG9nOxtHrqPA&e=.

-- Dr. Juliette Mammei

Associate Professor Dep't of Physics and Astronomy University of Manitoba Winnipeg, MB R3T 2N2 Canada

hanjie1 commented 3 years ago

Line 152 and 201 are using different structure functions. One is for inelastic, the other is for quasi-elastic.

Ok. This actuality makes me wonder if the phase space generator is what we want. Here are my thoughts.

We detect particles at specific (E', theta). Since the detectors have resolutions, the cross section of this event should be integrated over \delta E' and \delta \theta:

xs=d\sigma/(d\Omega dE') \delta E' sin(theta)2pi(\delta \theta) (\delta E') and (\delta \theta) are the detector resolutions, which should be a constant. In my opinion, they don't need to be multiplied in the generator. Then "weight" is what we want (need to multiply by the luminosity later).

Honestly I don't quite understand what is the purpose of using "phaseSpaceFactor". But seems the "V" in "remollGenpInelastic" is what people want:

double efmax = mpbeamE/(mp + beamE(1.0-cos(th)));; double V = (fPh_max - fPh_min) (cos(fTh_min) - cos(fTh_max)) efmax;

Alos, maybe we should make the \theta sampling consistent between different generators. Looks like in "remollGenpElastic", it uses Tyler's way to generate \theta.

wdconinc commented 3 years ago

(\delta E') and (\delta \theta) are the detector resolutions, which should be a constant. In my opinion, they don't need to be multiplied in the generator. Then "weight" is what we want (need to multiply by the luminosity later).

The detector resolution (or, more accurately described as detector acceptance which includes aspects of resolution) should not be included at the generator level. What introduces this in results is the Monte Carlo integration technique: more events accepted means larger acceptance. This integral only works if you normalize correctly by including the phase space factor. This integral essentially becomes: \Int{accepted events} d\sigma/(d\Omega dE') d\Omega dE' = \Sum{accepted events} d\sigma/(d\Omega dE') \Delta \Omega \DeltaE', where the Deltas are the total phase space over number of simulated events (the total number of simulated events is normalized out elsewhere, not in each event generator).

There needs to be a normalization to the phase space in E' that is sampled. Whether the expression for efmax is correct, I don't know, but there should be something there for each randomly picked kinematic variable.

Also, maybe we should make the \theta sampling consistent between different generators

Yes and no. It would be great if everyone could use the same sampling. But some formulas for differential cross section that we use are written in terms of theta_COM (Moller) and others are written in terms of theta_lab (e.g. pion generator). This will already change the weighting. Rewriting the differential cross section formulas for the same kinematic variables is not always easy (or possible, when we use models based on lab data or other external input, e.g. Konrad Aniol's hyperon generator). By 'normalizing' with the phase space factors into an event weight, we avoid the need for having to pick a single sampling strategy.

Another reason to avoid the same sampling is that some processes would be very inefficient to simulate that way since the cross section changes rapidly (at least in theory, I don't think we are actually affected by this unless maybe in the very forward direction). If we wanted to be efficient in simulation we would pick our sampling much more judiciously so as to reduce variance in the detectors we actually care about...

wdconinc commented 3 years ago

@cipriangal @hanjie1 Was this resolved in today's meeting? Is there a remaining issue here?

hanjie1 commented 3 years ago

Yes. It’s resolved.

On Thu, Mar 11, 2021 at 7:54 PM Wouter Deconinck @.***> wrote:

@cipriangal https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_cipriangal&d=DwMCaQ&c=CJqEzB1piLOyyvZjb8YUQw&r=NOYpK-Jh_cR_vuzyWVop-w&m=xrXMg3mdc47X0j7_cVBJB-3k6dMQ_1TdzyvuRwueqeI&s=wv-kTsD7KbHj1IMRPlVRlF_B4FGiHr5IbeNj_QCS_XY&e= @hanjie1 https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_hanjie1&d=DwMCaQ&c=CJqEzB1piLOyyvZjb8YUQw&r=NOYpK-Jh_cR_vuzyWVop-w&m=xrXMg3mdc47X0j7_cVBJB-3k6dMQ_1TdzyvuRwueqeI&s=1fsFV_UOWOAW_JA72IXivYSrYIdRHaQ9gVB7UW1TIrU&e= Was this resolved in today's meeting? Is there a remaining issue here?

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_JeffersonLab_remoll_issues_443-23issuecomment-2D797156786&d=DwMCaQ&c=CJqEzB1piLOyyvZjb8YUQw&r=NOYpK-Jh_cR_vuzyWVop-w&m=xrXMg3mdc47X0j7_cVBJB-3k6dMQ_1TdzyvuRwueqeI&s=4vrjpJrO8B8MGz5rNKdofRUu1Z7IJFJGbKWiCBpZxwc&e=, or unsubscribe https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AFSR25NGWHAD3EQG7YAW4U3TDFJ3TANCNFSM4YTWQS4A&d=DwMCaQ&c=CJqEzB1piLOyyvZjb8YUQw&r=NOYpK-Jh_cR_vuzyWVop-w&m=xrXMg3mdc47X0j7_cVBJB-3k6dMQ_1TdzyvuRwueqeI&s=SxwE-TYZIK3tp3UebBbYzZTP_4P_m4VXUH4SZKnKl5g&e= .