Only put particles which are actually transported on the stack.
[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"
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//*----------------------------------------------------------------------*
e3926667 34 static Double_t emassmev = PAPROP.am[9] * 1000.;
35
f827183d 36 Int_t lbhabh = 0;
37 if (EVTFLG.ldltry == 1) {
7dac99f1 38 if (EMFSTK.ichemf[EMFSTK.npemf-1] * EMFSTK.ichemf[EMFSTK.npemf-2] < 0) lbhabh = 1;
f827183d 39 }
e3926667 40
41// mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR
42// mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR
f827183d 43// EMFSTK.espark = spare real variables available for
44// EMFSTK.iespak = spare integer variables available for
e3926667 45// TRACKR.spausr = user defined spare variables for the current particle
46// TRACKR.ispusr = user defined spare flags for the current particle
47// EMFSTK.louemf = user flag
48// TRACKR.llouse = user defined flag for the current particle
f827183d 49
50 Int_t npnw, ispr;
e3926667 51 for (npnw = EMFSTK.npstrt-1; npnw <= EMFSTK.npemf-1; npnw++) {
f827183d 52
e3926667 53 for (ispr = 0; ispr <= mkbmx1-1; ispr++)
f827183d 54 EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr];
e3926667 55
56 for (ispr = 0; ispr <= mkbmx2-1; ispr++)
f827183d 57 EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr];
e3926667 58
f827183d 59 EMFSTK.louemf[npnw] = TRACKR.llouse;
60 }
61
62// Get the pointer to the VMC
fbf08100 63 TFluka* fluka = (TFluka*) gMC;
2bc4c610 64 Int_t verbosityLevel = fluka->GetVerbosityLevel();
65 Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
fbf08100 66 fluka->SetTrackIsNew(kTRUE);
e3926667 67
f827183d 68// Get the stack produced from the generator
69 TVirtualMCStack* cppstack = fluka->GetStack();
70
71// EVTFLG.ntrcks = track number
72// Increment the track number and put it into the last flag
73
74 Int_t kp;
e3926667 75
7dac99f1 76 for (kp = EMFSTK.npstrt - 1; kp <= EMFSTK.npemf - 1; kp++) {
e3926667 77
78// Ckeck transport cut first
79 Int_t ireg = EMFSTK.iremf[kp];
80 Double_t cut = (TMath::Abs(EMFSTK.ichemf[kp]) == 1) ? EFMRGN.ecut[ireg-1] : EFMRGN.pcut[ireg-1];
81 Double_t e = EMFSTK.etemf[kp];
82 if ((e < cut)
83 && (
84 (EMFSTK.ichemf[kp] == 0) ||
85 (EMFSTK.ichemf[kp] == -1) ||
86 (EMFSTK.ichemf[kp] == 1 && EFMRGN.pcut[ireg-1] > emassmev)
87 )
88 )
89 {
90 EMFSTK.iespak[kp][mkbmx2-1] = -1;
91 EMFSTK.iespak[kp][mkbmx2-2] = 0;
92 continue;
93 }
f827183d 94
95//* save the parent track number and reset it at each loop
610cd189 96 Int_t done = 0;
81ee7392 97 Int_t parent = TRACKR.ispusr[mkbmx2-1];
f827183d 98 Int_t flukaid = 0;
3b8c325d 99
e3926667 100// Identify particle type
101 if (EMFSTK.ichemf[kp] == -1) flukaid = 3;
102 else if (EMFSTK.ichemf[kp] == 0) flukaid = 7;
103 else if (EMFSTK.ichemf[kp] == 1) flukaid = 4;
104
105
106 e *= emvgev;
107 Int_t pdg = fluka->PDGFromId(flukaid);
d15554e0 108 Double_t p = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
109 Double_t px = p * EMFSTK.u[kp];
110 Double_t pz = p * EMFSTK.v[kp];
111 Double_t py = p * EMFSTK.w[kp];
112 Double_t tof = EMFSTK.agemf[kp];
113 Double_t polx = EMFSTK.upol[kp];
114 Double_t poly = EMFSTK.vpol[kp];
115 Double_t polz = EMFSTK.wpol[kp];
116 Double_t vx = EMFSTK.x[kp];
117 Double_t vy = EMFSTK.y[kp];
118 Double_t vz = EMFSTK.z[kp];
7dac99f1 119 Double_t weight = EMFSTK.wtemf[kp];
120
f827183d 121 Int_t ntr;
122 TMCProcess mech;
123 Int_t is = 0;
124
125//* case of no parent left (pair, photoelectric, annihilation):
126//* all secondaries are true
127 if ((EVTFLG.lpairp == 1) || (EVTFLG.lphoel == 1) ||
128 (EVTFLG.lannfl == 1) || (EVTFLG.lannrs == 1)) {
81ee7392 129
130 if (EVTFLG.lpairp == 1) mech = kPPair;
131 else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
132 else mech = kPAnnihilation;
642f15cf 133 cppstack->PushTrack(done, parent, pdg,
81ee7392 134 px, py, pz, e, vx, vy, vz, tof,
135 polx, poly, polz, mech, ntr, weight, is);
e3926667 136 if (debug) cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << "energy " << e-PAPROP.am[flukaid+6] << endl;
d15554e0 137
81ee7392 138 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
48aa0331 139 EMFSTK.iespak[kp][mkbmx2-2] = 0;
f827183d 140 } // end of lpairp, lphoel, lannfl, lannrs
141
142//* Compton: secondary is true only if charged (e+, e-)
143 else if ((EVTFLG.lcmptn == 1)) {
7dac99f1 144
145 if (EMFSTK.ichemf[kp] != 0) {
81ee7392 146 mech = kPCompton;
642f15cf 147 cppstack->PushTrack(done, parent, pdg,
81ee7392 148 px, py, pz, e, vx, vy, vz, tof,
149 polx, poly, polz, mech, ntr, weight, is);
2bc4c610 150 if (debug) cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
81ee7392 151 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
48aa0331 152 EMFSTK.iespak[kp][mkbmx2-2] = 0;
81ee7392 153 }
f827183d 154 } // end of lcmptn
81ee7392 155
f827183d 156//* Bremsstrahlung: true secondary only if charge = 0 (photon)
157 else if ((EVTFLG.lbrmsp == 1)) {
7dac99f1 158 if (EMFSTK.ichemf[kp] == 0) {
81ee7392 159 mech = kPBrem;
642f15cf 160 cppstack->PushTrack(done, parent, pdg,
81ee7392 161 px, py, pz, e, vx, vy, vz, tof,
162 polx, poly, polz, mech, ntr, weight, is);
2bc4c610 163 if (debug) cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
81ee7392 164 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
48aa0331 165 EMFSTK.iespak[kp][mkbmx2-2] = 0;
81ee7392 166 }
f827183d 167 } // end of lbrmsp
81ee7392 168
f827183d 169//* Delta ray: If Bhabha, true secondary only if negative (electron)
170 else if ((EVTFLG.ldltry == 1)) {
81ee7392 171 if (lbhabh == 1) {
7dac99f1 172 if (EMFSTK.ichemf[kp] == -1) {
81ee7392 173 mech = kPDeltaRay;
642f15cf 174 cppstack->PushTrack(done, parent, pdg,
81ee7392 175 px, py, pz, e, vx, vy, vz, tof,
176 polx, poly, polz, mech, ntr, weight, is);
177 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
48aa0331 178 EMFSTK.iespak[kp][mkbmx2-2] = 0;
2bc4c610 179 if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
81ee7392 180 } // end of Bhabha
d15554e0 181 } // lbhabh == 1
81ee7392 182
f827183d 183//* Delta ray: Otherwise Moller: true secondary is the electron with
184//* lower energy, which has been put higher in the stack
7dac99f1 185 else if (kp == EMFSTK.npemf-1) {
81ee7392 186 mech = kPDeltaRay;
642f15cf 187 cppstack->PushTrack(done, parent, pdg,
81ee7392 188 px, py, pz, e, vx, vy, vz, tof,
189 polx, poly, polz, mech, ntr, weight, is);
2bc4c610 190 if (debug) cout << endl << " !!! stupre (Moller) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
81ee7392 191 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
48aa0331 192 EMFSTK.iespak[kp][mkbmx2-2] = 0;
81ee7392 193 } // end of Delta ray
f827183d 194 } // end of ldltry
81ee7392 195
f827183d 196 } // end of loop
f827183d 197} // end of stupre
198} // end of extern "C"
199