# Re: [fluka-discuss]: how to include the divergence in the source.f / photons

From: Francesc Salvat-Pujol <francesc.salvat.pujol_at_cern.ch>
Date: Thu, 25 Jul 2019 19:29:14 +0200

Hey Alex,

If I get you right, you want to sample beam particle initial directions
with a rectangular divergence (i.e. cosine uniformly sampled in a given
interval) around an arbitrary direction theta_f and phi_f in your
modified source.f.

On paper, what one would ordinarily do is first sample the initial
direction with the desired rectangular divergence around the z axis and,
then, apply a rotation onto the final direction. So you're almost there.

- I believe for a rectangular beam the divergence is to be understood as
a full width. Thus, to sample the (cosine of the) polar angle one
would rather do

CDH = COS ( HLFHLF * DIV )
COSTH = ONEONE - ( ONEONE - CDH ) * FLRNDM ( TTT )

alternatively

COSTH = CDH + ( ONEONE - CDH ) * FLRNDM ( TTT )

where DIV is the divergence in radians, and the azimuthal angle you
sample uniformly from 0 to 2*pi like you already do. By the way, the
DIVBM variable should still have the value you passed in WHAT(3) of
particular case.

- You then construct T{X,Y,Z}FLK and ensure it is unit normalized. If
you chose to do smth a la TZFLK(NPFLKA) =SQRT(...) keep an eye on its
sign (for big DIV).

- The resulting T{X,Y,Z}FLK should be fine, but they are still sampled
around the z axis. Let d_i=(tx,ty,tz). To get the "final" direction,
d_f=(tx',ty',tz'), one would apply on d_i a rotation of angle theta_f
around the y axis, R(theta_f,y), and, next, a rotation of angle phi_f
around the z axis, R(phi_f,z).

To be explicit, see the (hopefully not too far off) attached recipe,
without any particular effort to be numerically robust (which is a bad
idea... and exactly the reason why director cosines are preferrable
instead of angles...), but it will get you started. At the very least,
make sure that the final T{X,Y,Z}FLK components are normalized (such
that the sum of their squares is ONEONE).

Cesc

On Tue, Jul 23 2019, at 18:36 +0000, Ševčik Aleksandras wrote:
>
>Thank You, Sougata,
>
>I meant BEAMPOS WHAT(4) & (5) , i.e. the ones which gives the direction cosines of the beam itself.
>The divergence works itself great when beam is straight, with BEAMPOS WHAT(4)&(5) = 0...
>
>In case anyone can help with this question, also with the coding for the rectangular beam divergence, kindly please,
>
>Regards
>Alex
>
>-----Original Message-----
>From: Sougata <sougata_at_barc.gov.in>
>Sent: Tuesday, July 23, 2019 13:02
>To: Ševčik Aleksandras <aleksandras.sevcik_at_ktu.edu>
>Cc: fluka-discuss_at_fluka.org; owner-fluka-discuss_at_mi.infn.it
>Subject: RE: [fluka-discuss]: how to include the divergence in the source.f / photons
>
>Dear Alex,
>
>Any parameter from source card may be passed to source.f file through whasou(#).
>
>e.g
>Coding for 200 mrad divergence beam where WHAT(3)=0.98007 im source card.
>
> * Cosines (tx,ty,tz)
> COSTH = ONEONE-(1-whasou(3))*FLRNDM(TTT)
> IF(COSTH.GT.ONEONE) COSTH=ONEONE
> IF(COSTH.LT.-ONEONE) COSTH=-ONEONE
> PHI = TWOPIP*FLRNDM(SSS)
> UBEAM=SQRT(ONEONE-COSTH*COSTH)*COS(PHI)
> VBEAM=SQRT(ONEONE-COSTH*COSTH)*SIN(PHI)
> WBEAM=COSTH
> TXFLK (NPFLKA) = UBEAM
> TYFLK (NPFLKA) = VBEAM
> * TZFLK (NPFLKA) = WBEAM
> TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
> & - TYFLK (NPFLKA)**2 )
>
>
>Regarding your query (2). I am not so experienced. Better consult some expert person.
>
>Regards,
>--
>SOUGATA RAKSHIT / सौगत रक्षित
>Scientific Officer(D) / वैज्ञानिक आधिकारीक(D) RSSD / आर एस एस डि Bhabha Atomic Research Centre / भाभा परमानु अनुसंधान केंद्र
>Mumbai-400 085 / मुंबई-400 085
>Tel: +91-022-2559 2214/ दूरभाष: +91-022-2559 2214
>
>On 2019-07-23 12:57, Ševčik Aleksandras wrote:
>> Dear Sougata Rakshit,
>>
>> Wonderful, it works indeed!
>> But may I ask two more moments here:
>>
>> 1) I need to set up the certain values ( 0.2 ) for BEAMPOS WHAT(4) &
>> (5) - direction cosine of the beam with respect to the x & y-axis. It
>> does not take values from the card now when I use this modified
>> routine. What additional line should be included for this?
>> 2) also this code works for angular flat beam only, but not for
>> rectangular flat beam - does the code differs much in such case?
>>
>> If you or anybody else could find an opportunity to help with this, I
>> will greatly appreciate such help,
>>
>> Regards
>> Alex
>>
>> -----Original Message-----
>> From: Sougata <sougata_at_barc.gov.in>
>> Sent: Tuesday, July 23, 2019 07:47
>> To: Ševčik Aleksandras <aleksandras.sevcik_at_ktu.edu>
>> Cc: fluka-discuss_at_fluka.org; owner-fluka-discuss_at_mi.infn.it
>> Subject: Re: [fluka-discuss]: how to include the divergence in the
>> source.f / photons
>>
>>
>> Dear Alex,
>>
>>
>> So we have to randomize direction cosine along beam axis from 0.98007
>> to 1 (i.e. cos (0 degree) to cos (200 mrad))
>>
>> you can try following code as
>>
>> * Cosines (tx,ty,tz)
>> COSTH = ONEONE-(1-0.98007)*FLRNDM(TTT)
>> IF(COSTH.GT.ONEONE) COSTH=ONEONE
>> IF(COSTH.LT.-ONEONE) COSTH=-ONEONE
>> PHI = TWOPIP*FLRNDM(SSS)
>> UBEAM=SQRT(ONEONE-COSTH*COSTH)*COS(PHI)
>> VBEAM=SQRT(ONEONE-COSTH*COSTH)*SIN(PHI)
>> WBEAM=COSTH
>> TXFLK (NPFLKA) = UBEAM
>> TYFLK (NPFLKA) = VBEAM
>> * TZFLK (NPFLKA) = WBEAM
>> TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
>> & - TYFLK (NPFLKA)**2 )
>>
>> I think it will work.
>>
>> Regards,
>> --
>> SOUGATA RAKSHIT / सौगत रक्षित
>> Scientific Officer(D) / वैज्ञानिक आधिकारीक(D) RSSD / आर एस एस डि
>> Bhabha Atomic Research Centre / भाभा परमानु अनुसंधान केंद्र
>> Mumbai-400 085 / मुंबई-400 085
>> Tel: +91-022-2559 2214/ दूरभाष: +91-022-2559 2214
>>
>>
>> On 2019-07-23 00:21, Ševčik Aleksandras wrote:
>>> Dear experts,
>>>
>>> I am using the attached source.f to sample the designed photon
>>> spectrum, however, I noticed that BEAM WHAT(3) does not respond to
>>> the values if source routine is activated - just as stated in the manual.
>>>
>>> Is there any "simple" way to make my source beam with divergence 200
>>> mrad, or only additional piece of code must be included in the
>>> source.f?
>>>
>>> Looking through old forum archives I found similar code for that
>>> purpose:
>>>
>>> Divbm_Sigma = Divbm*0.2
>>>
>>> CALL FLNRR2 (RGAUSSDIVX,RGAUSSDIVY)
>>>
>>> DIV_VECTORX = RGAUSSDIVX*Divbm_Sigma
>>>
>>> DIV_VECTORY = RGAUSSDIVY*Divbm_Sigma
>>>
>>> TXHLP = TAN( DIV_VECTORX )
>>>
>>> TYHLP = TAN( DIV_VECTORY )
>>>
>>> THELP = SQRT( TXHLP*TXHLP + TYHLP*TYHLP + ONEONE )
>>>
>>> TXFLK (NPFLKA) = TXHLP / THELP
>>>
>>> TYFLK (NPFLKA) = TYHLP / THELP
>>>
>>> TZFLK (NPFLKA) = SQRT ( ONEONE - LTXFLK (NPFLKA)**2 - TYFLK
>>> (NPFLKA)**2 )
>>>
>>> TXFLK (NPFLKA) = UBEAM
>>>
>>> TYFLK (NPFLKA) = VBEAM
>>>
>>> TZFLK (NPFLKA) = WBEAM
>>>
>>> TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
>>>
>>> & - TYFLK (NPFLKA)**2 )
>>>
>>> But simply applying it the code is not compiled, and my lack of
>>> experience with user routines is a big obstacle to move forward from
>>> this point. In addition, my BEAMPOS is defined with WHAT(1)-(WHAT(3)
>>> and also rotation value for WHAT(5) is used. Also the beam is
>>> directed to positive Z direction.
>>>
>>> Could anyone share a similar working source.f file which would
>>> directly use BEAM WhHAT(3) value, or maybe guide me through this
>>> coding problem?
>>>
>>> Both source.f routines - the original and the one I tried to modify -
>>> are attached.
>>>
>>> Thank you very much,
>>>
>>> Rgds
>>>
>>> Alex
>>
>>
>>
>> ______________________________________________________________________
>> ____ You can manage unsubscription from this mailing list at
>> https://www.fluka.org/fluka.php?id=acc_info
>
>
>
>__________________________________________________________________________
>You can manage unsubscription from this mailing list at https://www.fluka.org/fluka.php?id=acc_info
>

```--
Francesc Salvat Pujol
CERN-EN/STI
CH-1211 Geneva 23
Switzerland
Tel: +41 22 76 64011
Fax: +41 22 76 69474

```

__________________________________________________________________________
You can manage unsubscription from this mailing list at https://www.fluka.org/fluka.php?id=acc_info

Received on Thu Jul 25 2019 - 23:07:18 CEST

This archive was generated by hypermail 2.3.0 : Thu Jul 25 2019 - 23:07:31 CEST