Correct energy of interacting particle after BREMS, DELTA, COMPTON
[u/mrichter/AliRoot.git] / TFluka / stupre.cxx
CommitLineData
f827183d 1#include <Riostream.h>
f827183d 2#ifndef WIN32
3# define stupre stupre_
4#else
5# define stupre STUPRE
6#endif
7//
8// Fluka include
9#include "Fdimpar.h" //(DIMPAR) fluka include
e3926667 10
f827183d 11// Fluka commons
12#include "Fdblprc.h" //(DBLPRC) fluka common
13#include "Femfstk.h" //(EMFSTK) fluka common
14#include "Fevtflg.h" //(EVTFLG) fluka common
15#include "Fpaprop.h" //(PAPROP) fluka common
16#include "Ftrackr.h" //(TRACKR) fluka common
e3926667 17#include "Femfrgn.h" //(EFMRGN) fluka common
f827183d 18
19//Virtual MC
20#include "TFluka.h"
d566901f 21#include "TFlukaCodes.h"
f827183d 22#include "TVirtualMCStack.h"
23#include "TVirtualMCApplication.h"
24#include "TParticle.h"
25#include "TVector3.h"
26
27extern "C" {
28void stupre()
29{
30//*----------------------------------------------------------------------*
31//* *
32//* SeT User PRoperties for Emf particles *
33//* *
34//*----------------------------------------------------------------------*
d566901f 35// Get the pointer to the VMC
36 TFluka* fluka = (TFluka*) gMC;
37 static Double_t emassmev = PAPROP.am[kFLUKAelectron + 6] * 1000.;
e3926667 38
f827183d 39 Int_t lbhabh = 0;
40 if (EVTFLG.ldltry == 1) {
7dac99f1 41 if (EMFSTK.ichemf[EMFSTK.npemf-1] * EMFSTK.ichemf[EMFSTK.npemf-2] < 0) lbhabh = 1;
f827183d 42 }
d566901f 43
e3926667 44// mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR
45// mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR
f827183d 46// EMFSTK.espark = spare real variables available for
47// EMFSTK.iespak = spare integer variables available for
e3926667 48// TRACKR.spausr = user defined spare variables for the current particle
49// TRACKR.ispusr = user defined spare flags for the current particle
50// EMFSTK.louemf = user flag
51// TRACKR.llouse = user defined flag for the current particle
f827183d 52
53 Int_t npnw, ispr;
e3926667 54 for (npnw = EMFSTK.npstrt-1; npnw <= EMFSTK.npemf-1; npnw++) {
4aba9d66 55
d566901f 56 for (ispr = 0; ispr <= mkbmx1-1; ispr++)
4aba9d66 57 EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr];
58
d566901f 59 for (ispr = 0; ispr <= mkbmx2-1; ispr++)
4aba9d66 60 EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr];
61
d566901f 62 EMFSTK.louemf[npnw] = TRACKR.llouse;
f827183d 63 }
64
d566901f 65
2bc4c610 66 Int_t verbosityLevel = fluka->GetVerbosityLevel();
67 Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
82a3f706 68 debug = 100;
69
fbf08100 70 fluka->SetTrackIsNew(kTRUE);
e3926667 71
f827183d 72// Get the stack produced from the generator
73 TVirtualMCStack* cppstack = fluka->GetStack();
74
75// EVTFLG.ntrcks = track number
76// Increment the track number and put it into the last flag
77
78 Int_t kp;
e3926667 79
7dac99f1 80 for (kp = EMFSTK.npstrt - 1; kp <= EMFSTK.npemf - 1; kp++) {
e3926667 81
4aba9d66 82// Check transport cut first
d566901f 83 Int_t ireg = EMFSTK.iremf[kp];
84 Double_t cut = (TMath::Abs(EMFSTK.ichemf[kp]) == 1) ? EMFRGN.elethr[ireg-1] : EMFRGN.phothr[ireg-1];
e3926667 85 Double_t e = EMFSTK.etemf[kp];
b7940702 86
e3926667 87 if ((e < cut)
4aba9d66 88 && (
89 (EMFSTK.ichemf[kp] == 0) ||
90 (EMFSTK.ichemf[kp] == -1) ||
91 (EMFSTK.ichemf[kp] == 1 && EMFRGN.phothr[ireg-1] > emassmev)
92 )
93 )
e3926667 94 {
4aba9d66 95 EMFSTK.iespak[kp][mkbmx2-1] = TRACKR.ispusr[mkbmx2-1];
96 EMFSTK.iespak[kp][mkbmx2-2] = 0;
97 EMFSTK.iespak[kp][mkbmx2-3] = TRACKR.jtrack;
98 EMFSTK.irlatt[kp] = TRACKR.lt1trk;
99 continue;
e3926667 100 }
f827183d 101
d566901f 102// Save the parent track number and reset it at each loop
610cd189 103 Int_t done = 0;
81ee7392 104 Int_t parent = TRACKR.ispusr[mkbmx2-1];
f827183d 105 Int_t flukaid = 0;
4aba9d66 106
e3926667 107// Identify particle type
d566901f 108 if (EMFSTK.ichemf[kp] == -1) flukaid = kFLUKAelectron;
109 else if (EMFSTK.ichemf[kp] == 0) flukaid = kFLUKAphoton;
110 else if (EMFSTK.ichemf[kp] == 1) flukaid = kFLUKApositron;
4aba9d66 111
e3926667 112
d566901f 113 e *= emvgev;
e3926667 114 Int_t pdg = fluka->PDGFromId(flukaid);
d15554e0 115 Double_t p = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
81f1d030 116 Double_t px = p * EMFSTK.uemf[kp];
b7940702 117 Double_t py = p * EMFSTK.vemf[kp];
118 Double_t pz = p * EMFSTK.wemf[kp];
d15554e0 119 Double_t tof = EMFSTK.agemf[kp];
120 Double_t polx = EMFSTK.upol[kp];
121 Double_t poly = EMFSTK.vpol[kp];
122 Double_t polz = EMFSTK.wpol[kp];
81f1d030 123 Double_t vx = EMFSTK.xemf[kp];
124 Double_t vy = EMFSTK.yemf[kp];
125 Double_t vz = EMFSTK.zemf[kp];
7dac99f1 126 Double_t weight = EMFSTK.wtemf[kp];
127
f827183d 128 Int_t ntr;
129 TMCProcess mech;
130 Int_t is = 0;
131
132//* case of no parent left (pair, photoelectric, annihilation):
133//* all secondaries are true
134 if ((EVTFLG.lpairp == 1) || (EVTFLG.lphoel == 1) ||
135 (EVTFLG.lannfl == 1) || (EVTFLG.lannrs == 1)) {
4aba9d66 136
137 if (EVTFLG.lpairp == 1) mech = kPPair;
138 else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
139 else mech = kPAnnihilation;
642f15cf 140 cppstack->PushTrack(done, parent, pdg,
4aba9d66 141 px, py, pz, e, vx, vy, vz, tof,
142 polx, poly, polz, mech, ntr, weight, is);
143 if (debug) cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << " energy " << e-PAPROP.am[flukaid+6] << endl;
d15554e0 144
4aba9d66 145 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
146 EMFSTK.iespak[kp][mkbmx2-2] = 0;
f827183d 147 } // end of lpairp, lphoel, lannfl, lannrs
4aba9d66 148
f827183d 149//* Compton: secondary is true only if charged (e+, e-)
150 else if ((EVTFLG.lcmptn == 1)) {
7dac99f1 151
4aba9d66 152 if (EMFSTK.ichemf[kp] != 0) {
153 mech = kPCompton;
154 cppstack->PushTrack(done, parent, pdg,
155 px, py, pz, e, vx, vy, vz, tof,
156 polx, poly, polz, mech, ntr, weight, is);
157 if (debug) cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
158 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
159 EMFSTK.iespak[kp][mkbmx2-2] = 0;
82a3f706 160 } else {
161 fluka->SetPint(px, py, pz, e);
162 }
f827183d 163 } // end of lcmptn
4aba9d66 164
f827183d 165//* Bremsstrahlung: true secondary only if charge = 0 (photon)
166 else if ((EVTFLG.lbrmsp == 1)) {
4aba9d66 167 if (EMFSTK.ichemf[kp] == 0) {
168 mech = kPBrem;
169 cppstack->PushTrack(done, parent, pdg,
170 px, py, pz, e, vx, vy, vz, tof,
171 polx, poly, polz, mech, ntr, weight, is);
172 if (debug) cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
173 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
174 EMFSTK.iespak[kp][mkbmx2-2] = 0;
82a3f706 175 } else {
176 fluka->SetPint(px, py, pz, e);
177 }
178
f827183d 179 } // end of lbrmsp
4aba9d66 180
f827183d 181//* Delta ray: If Bhabha, true secondary only if negative (electron)
182 else if ((EVTFLG.ldltry == 1)) {
4aba9d66 183 if (lbhabh == 1) {
184 if (EMFSTK.ichemf[kp] == -1) {
185 mech = kPDeltaRay;
186 cppstack->PushTrack(done, parent, pdg,
82a3f706 187 px, py, pz, e, vx, vy, vz, tof,
188 polx, poly, polz, mech, ntr, weight, is);
4aba9d66 189 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
190 EMFSTK.iespak[kp][mkbmx2-2] = 0;
82a3f706 191 if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
192 } else {
193 fluka->SetPint(px, py, pz, e);
194 }
4aba9d66 195 } // lbhabh == 1
196
f827183d 197//* Delta ray: Otherwise Moller: true secondary is the electron with
198//* lower energy, which has been put higher in the stack
4aba9d66 199 else if (kp == EMFSTK.npemf-1) {
200 mech = kPDeltaRay;
201 cppstack->PushTrack(done, parent, pdg,
202 px, py, pz, e, vx, vy, vz, tof,
203 polx, poly, polz, mech, ntr, weight, is);
204 if (debug) cout << endl << " !!! stupre (Moller) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
205 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
206 EMFSTK.iespak[kp][mkbmx2-2] = 0;
82a3f706 207 } else {
208 fluka->SetPint(px, py, pz, e);
209 }
f827183d 210 } // end of ldltry
4aba9d66 211
f827183d 212 } // end of loop
f827183d 213} // end of stupre
214} // end of extern "C"
215