Energy scoring along step in user routine, biasing, and particle weight

Advanced Features and User Routines
1

Dear colleagues,
somehow following up on the question that was discussed here, I am now seeing an unexpected (for me) effect related to energy scoring along a particle step when quenching is activated by means of the USERDUMP card, SDUM=UDQUENCH.

Specifically, in my setup I am scoring energy deposition by muons passing through a 4-cm thick plastic scintillator, with first Birks parameter set to 0.013. I do so via mgdraw.f

If I set the weight of the primary muon to one (via a custom source.f code), I see a proper energy deposition spectrum, that is almost the same as that I obtain if I de-activate the quenching mechanism (as expected).

However, If I set the weight of the primary to 1E-3, I see the following:

  • If quenching is not activated, the energy spectrum looks as the one with primary weight set to 1.
  • If quenching is activated, instead, all events have a very small energy deposition.

Is there a correlation between the deposited energy (with quenching) and the weight of the particle?

The code I am using in mgdraw.f to handle quenching is:

*  +-------------------------------------------------------------------*
*  |  Quenching is activated
      IF ( LQEMGD ) THEN
         IF ( MTRACK .GT. 0 ) THEN
           RULLL  = ZERZER
           CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN )
         END IF
      END IF
*  |  End of  quenching
*  +-------------------------------------------------------------------* 
	    
	      CALL digicontinuos(JTRACK,ETRACK,MTRACK,DTRACK,DTQUEN,
     &	      XTRACK,YTRACK,ZTRACK,MREG,LT1TRK,LT2TRK,MNAME)
	  end if

digicontinuos is the function I use to score at the event level the energy deposition. As I reported before, in case the primary weight of the particle is set to one, everything works as expected.

Thanks,
Andrea

2

Hi Andrea,
as a general word of caution, event-by-event scoring - as required to get an energy deposition spectrum - is somewhat incompatible with biasing, since the latter breaks natural correlations.
Said so, the wrong outcome you usefully report should be fixed by adding the following loops before and after the QUENMG call:

            DO I = 1, MTRACK
               DTRACK (I) = DTRACK (I) * WSCRNG
            END DO
            CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN )
            DO I = 1, MTRACK
               DTRACK (I) = DTRACK (I) / WSCRNG
               DO JBK = 1, NQEMGD
                  DTQUEN (I,JBK) = DTQUEN (I,JBK) / WSCRNG
               END DO
            END DO

Let us know in case of persistent issue.

3

Dear Francesco,
thanks for your response. I tried to implement this fix and it works, thanks!!

I agree that, in general, it is very dangerous to use biasing and then score observables on an event-by-event basis, since biasing breaks natural correlations.

In my case, I think this is under control. I am interested in the production of mu+ mu- pairs from gamma conversion in a target where a 100 GeV e- impinges. The e- starts an EM shower, and one of the photons of the shower can generate a mu+ mu- pair.

The bias factor that I am using for this process (i.e. the factor I increase the cross section for) is low enough to keep this as a “rare” process. Furthermore, following the suggestion I got in another topic, I am killing the initial state photon.

Thanks again,
Bests,
Andrea

1 Like
6

Hi,
I was too fast in replying that the proposed solution was working correctly, since I did just a quick check looking at few events. The plot below, however, shows you the energy distribution of 100 GeV/c muons in a 5-cm thick plastic scintillator with quenching activated via USERDUMP, SDUM=UDQUENCH (quenching factor 0.0126 g cm2 / MeV).

In blue, the result if I set, as usual, the weight of the primary muon to 1.

In red, the result if I set this weight to 0.002 (via source.f): the weight is the same for all primaries.

As you see, now the order of magnitude is the same, but the distribution looks significantly different. By the way, I also tried to implement “manually” quenching in mgdraw.f, by the formula E’ = E * 1 / (1+kB * dE/dX), where “dE/dX” is computed taking the energy deposition along the step divided by the step length. The result is the same as the red one below.

Thanks,

Andrea

image
image1280×776 98.2 KB

1 Like
7

I understand that, to build the above spectra, you count any energy deposition occurrence as 1, rather than WSCRNG, but you expect the latter to be always 0.002 in case of biasing, being then just a scaling factor that does not alter the red spectrum shape.
For convenience, can you please upload or make available your input and the source, mgdraw and digicontinuos routines?

8

Dear Francesco,
this is correct - since the weight is always 0.002, then the result should be the same as if the weight is one for all events.

I attach here the codes that you asked for: the input, and all the other user files. To compile the executable rootfluka, you can use the makefile in the tar file.

The output of the simulation is a text file (CaloHits.d), reporting the energy deposited per each event.
The relevant lines to look for are:

EVENT

this marks the beginning of an event

VETO
E1 E2 E3

where “E1”, “E2”, “E3” are three doubles. E2 is the energy deposited in the 5-cm thick counter.
Bests,
Andrea

fluka.tgz (53.3 KB)

www.ge.infn.it/~celentan/flukaCodes/flukaCodes.tar

9

Dear Francesco,
I removed the (wrong) corrections to the ENDRAW routine, but still I see a difference when I set the weight of primary particles to 0.002.

10

Thanks Andrea, and sorry for the slow convergence.
My lines above were indeed meant for the first QUENMG call only, dealing with energy depositions along the ionization track. But point-like energy depositions handled in the ENDRAW entry also required a (different) fix, namely:

         CALL QUENMG ( ICODE, MREG, RULLL, DTQUEN )
         DO JBK = 1, NQEMGD
            DTQUEN (1,JBK) = DTQUEN (1,JBK) / WSCRNG
         END DO
         RULLL = RULLL / WSCRNG

[mysterious red color not intended].
Moreover, in your digisingle routine, you should make use of RULLL, and not RULL.
All these additional lines will no longer be needed from the next release, as we’ll cure the issue directly inside the QUENMG code.
Of course, the last word remains yours, to make sure that no enduring problem is left.

11

Dear Francesco,
thanks. I implemented this new fix. The result is reported below (as before: blue means with weight=1 for all primaries, red means weight=0.002 for all primaries). It looks now everything works!

image

For my curiosity, may I ask you how the particle weight enters the routine for the quenching in QUENMG?

Thanks,
Bests,
Andrea

12

That’s exactly the point. For the purposes of event-by-event energy deposition scoring, the statistical weight should not be embedded in the deposited quantity (horizontal axis of your spectrum), rather separately kept for occurrence counting (vertical axis) - being always careful about the applicability of event-by-event scoring in the presence of biasing. In the case of energy deposition along the track, the relevant energy array does not embed the statistical weight, but the quenmg routine erroneously meant to remove it (you can interpret my temporary fix in mgdraw accordingly). In the case of point-like energy deposition, the RULL quantity does contain the statistical weight, which quenmg has to appropriately remove (it will eventually do it and return back a clean RULLL quantity, as well as DTQUEN).