Parent id and current track id handling corrected.
[u/mrichter/AliRoot.git] / TFluka / stupre.cxx
CommitLineData
f827183d 1#include <Riostream.h>
2#include "AliRun.h"
3#include "TFluka.h"
4#ifndef WIN32
5# define stupre stupre_
6#else
7# define stupre STUPRE
8#endif
9//
10// Fluka include
11#include "Fdimpar.h" //(DIMPAR) fluka include
12// Fluka commons
13#include "Fdblprc.h" //(DBLPRC) fluka common
14#include "Femfstk.h" //(EMFSTK) fluka common
15#include "Fevtflg.h" //(EVTFLG) fluka common
16#include "Fpaprop.h" //(PAPROP) fluka common
17#include "Ftrackr.h" //(TRACKR) fluka common
18
19//Virtual MC
20#include "TFluka.h"
21#include "TVirtualMCStack.h"
22#include "TVirtualMCApplication.h"
23#include "TParticle.h"
24#include "TVector3.h"
25
26extern "C" {
27void stupre()
28{
29//*----------------------------------------------------------------------*
30//* *
31//* SeT User PRoperties for Emf particles *
32//* *
33//*----------------------------------------------------------------------*
34
35 Int_t lbhabh = 0;
36 if (EVTFLG.ldltry == 1) {
37 if (EMFSTK.iq[EMFSTK.np-1] * EMFSTK.iq[EMFSTK.np-2] < 0) lbhabh = 1;
38 }
39
40// mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR
41// mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR
42// EMFSTK.espark = spare real variables available for
43// EMFSTK.iespak = spare integer variables available for
44// TRACKR.spausr = user defined spare variables for the current particle
45// TRACKR.ispusr = user defined spare flags for the current particle
46// EMFSTK.louemf = user flag
47// TRACKR.llouse = user defined flag for the current particle
48
49 Int_t npnw, ispr;
50 for (npnw=EMFSTK.npstrt-1; npnw<=EMFSTK.np-1; npnw++) {
51
52 for (ispr=0; ispr<=mkbmx1-1; ispr++)
53 EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr];
54
55 for (ispr=0; ispr<=mkbmx2-1; ispr++)
56 EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr];
57
58 EMFSTK.louemf[npnw] = TRACKR.llouse;
59 }
60
61// Get the pointer to the VMC
62 TVirtualMC* fluka = TFluka::GetMC();
63// Get the stack produced from the generator
64 TVirtualMCStack* cppstack = fluka->GetStack();
65
66// EVTFLG.ntrcks = track number
67// Increment the track number and put it into the last flag
68
69 Int_t kp;
81ee7392 70 for (kp = EMFSTK.npstrt-1; npnw <= EMFSTK.np-1; kp++) {
f827183d 71
72//* save the parent track number and reset it at each loop
81ee7392 73 Int_t done = 1;
74
75 Int_t parent = TRACKR.ispusr[mkbmx2-1];
f827183d 76
f827183d 77 Int_t flukaid = 0;
78 if (EMFSTK.iq[kp] == -1) flukaid = 3;
79 else if (EMFSTK.iq[kp] == 0) flukaid = 7;
80 else if (EMFSTK.iq[kp] == 0) flukaid = 4;
81 Int_t pdg = fluka->PDGFromId(flukaid);
82 Double_t e = EMFSTK.e[kp]*emvgev;
83 Double_t p = sqrt(e*e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
84 Double_t px = p * EMFSTK.u[kp];
85 Double_t pz = p * EMFSTK.v[kp];
86 Double_t py = p * EMFSTK.w[kp];
87 Double_t tof = EMFSTK.agemf[kp];
88 Double_t polx = EMFSTK.upol[kp];
89 Double_t poly = EMFSTK.vpol[kp];
90 Double_t polz = EMFSTK.wpol[kp];
91 Double_t vx = EMFSTK.x[kp];
92 Double_t vy = EMFSTK.y[kp];
93 Double_t vz = EMFSTK.z[kp];
94 Double_t weight = EMFSTK.wt[kp];
95 Int_t ntr;
96 TMCProcess mech;
97 Int_t is = 0;
98
99//* case of no parent left (pair, photoelectric, annihilation):
100//* all secondaries are true
101 if ((EVTFLG.lpairp == 1) || (EVTFLG.lphoel == 1) ||
102 (EVTFLG.lannfl == 1) || (EVTFLG.lannrs == 1)) {
81ee7392 103
104 if (EVTFLG.lpairp == 1) mech = kPPair;
105 else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
106 else mech = kPAnnihilation;
f827183d 107 cppstack->SetTrack(done, parent, pdg,
81ee7392 108 px, py, pz, e, vx, vy, vz, tof,
109 polx, poly, polz, mech, ntr, weight, is);
110
111 cout << endl << " !!! stupre: ntr=" << ntr << " parent=" << parent << endl;
112 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
f827183d 113 } // end of lpairp, lphoel, lannfl, lannrs
114
115//* Compton: secondary is true only if charged (e+, e-)
116 else if ((EVTFLG.lcmptn == 1)) {
81ee7392 117 if (EMFSTK.iq[kp] != 0) {
118 mech = kPCompton;
119 cppstack->SetTrack(done, parent, pdg,
120 px, py, pz, e, vx, vy, vz, tof,
121 polx, poly, polz, mech, ntr, weight, is);
122 cout << endl << " !!! stupre: ntr=" << ntr << " parent=" << parent << endl;
123 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
124 }
f827183d 125 } // end of lcmptn
81ee7392 126
f827183d 127//* Bremsstrahlung: true secondary only if charge = 0 (photon)
128 else if ((EVTFLG.lbrmsp == 1)) {
81ee7392 129 if (EMFSTK.iq[kp] == 0) {
130 mech = kPBrem;
131 cppstack->SetTrack(done, parent, pdg,
132 px, py, pz, e, vx, vy, vz, tof,
133 polx, poly, polz, mech, ntr, weight, is);
134 cout << endl << " !!! stupre: ntr=" << ntr << " parent=" << parent << endl;
135 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
136 }
f827183d 137 } // end of lbrmsp
81ee7392 138
f827183d 139//* Delta ray: If Bhabha, true secondary only if negative (electron)
140 else if ((EVTFLG.ldltry == 1)) {
81ee7392 141 if (lbhabh == 1) {
142 if (EMFSTK.iq[kp] == 0) {
143 mech = kPDeltaRay;
144 cppstack->SetTrack(done, parent, pdg,
145 px, py, pz, e, vx, vy, vz, tof,
146 polx, poly, polz, mech, ntr, weight, is);
147 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
148 } // end of Bhabha
149 }
150
f827183d 151//* Delta ray: Otherwise Moller: true secondary is the electron with
152//* lower energy, which has been put higher in the stack
81ee7392 153 else if (kp == EMFSTK.np) {
154 mech = kPDeltaRay;
155 cppstack->SetTrack(done, parent, pdg,
156 px, py, pz, e, vx, vy, vz, tof,
157 polx, poly, polz, mech, ntr, weight, is);
158 cout << endl << " !!! stupre: ntr=" << ntr << " parent=" << parent << endl;
159 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
160 } // end of Delta ray
f827183d 161 } // end of ldltry
81ee7392 162
f827183d 163 } // end of loop
81ee7392 164
f827183d 165// !!! TO BE CONFIRMED !!!
166 EVTFLG.ntrcks = EMFSTK.iespak[EMFSTK.np-1][mkbmx2-1];
81ee7392 167
f827183d 168} // end of stupre
169} // end of extern "C"
170