3 # define stupre stupre_
9 #include "Fdimpar.h" //(DIMPAR) fluka include
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
17 #include "Femfrgn.h" //(EFMRGN) fluka common
21 #include "TFlukaCodes.h"
22 #include "TVirtualMCStack.h"
23 #include "TVirtualMCApplication.h"
24 #include "TParticle.h"
30 //*----------------------------------------------------------------------*
32 //* SeT User PRoperties for Emf particles *
34 //*----------------------------------------------------------------------*
35 // Get the pointer to the VMC
36 TFluka* fluka = (TFluka*) gMC;
37 static Double_t emassmev = PAPROP.am[kFLUKAelectron + 6] * 1000.;
40 if (EVTFLG.ldltry == 1) {
41 if (EMFSTK.ichemf[EMFSTK.npemf-1] * EMFSTK.ichemf[EMFSTK.npemf-2] < 0) lbhabh = 1;
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
46 // EMFSTK.espark = spare real variables available for
47 // EMFSTK.iespak = spare integer variables available for
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
54 for (npnw = EMFSTK.npstrt-1; npnw <= EMFSTK.npemf-1; npnw++) {
56 for (ispr = 0; ispr <= mkbmx1-1; ispr++)
57 EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr];
59 for (ispr = 0; ispr <= mkbmx2-1; ispr++)
60 EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr];
62 EMFSTK.louemf[npnw] = TRACKR.llouse;
66 Int_t verbosityLevel = fluka->GetVerbosityLevel();
67 Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
69 fluka->SetTrackIsNew(kTRUE);
71 // Get the stack produced from the generator
72 TVirtualMCStack* cppstack = fluka->GetStack();
74 // EVTFLG.ntrcks = track number
75 // Increment the track number and put it into the last flag
79 for (kp = EMFSTK.npstrt - 1; kp <= EMFSTK.npemf - 1; kp++) {
81 // Check transport cut first
82 Int_t ireg = EMFSTK.iremf[kp];
83 Double_t cut = (TMath::Abs(EMFSTK.ichemf[kp]) == 1) ? EMFRGN.elethr[ireg-1] : EMFRGN.phothr[ireg-1];
84 Double_t e = EMFSTK.etemf[kp];
88 (EMFSTK.ichemf[kp] == 0) ||
89 (EMFSTK.ichemf[kp] == -1) ||
90 (EMFSTK.ichemf[kp] == 1 && EMFRGN.phothr[ireg-1] > emassmev)
94 EMFSTK.iespak[kp][mkbmx2-1] = TRACKR.ispusr[mkbmx2-1];
95 EMFSTK.iespak[kp][mkbmx2-2] = 0;
96 EMFSTK.iespak[kp][mkbmx2-3] = TRACKR.jtrack;
97 EMFSTK.irlatt[kp] = TRACKR.lt1trk;
101 // Save the parent track number and reset it at each loop
103 Int_t parent = TRACKR.ispusr[mkbmx2-1];
106 // Identify particle type
107 if (EMFSTK.ichemf[kp] == -1) flukaid = kFLUKAelectron;
108 else if (EMFSTK.ichemf[kp] == 0) flukaid = kFLUKAphoton;
109 else if (EMFSTK.ichemf[kp] == 1) flukaid = kFLUKApositron;
113 Int_t pdg = fluka->PDGFromId(flukaid);
114 Double_t p = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
115 Double_t px = p * EMFSTK.uemf[kp];
116 Double_t py = p * EMFSTK.vemf[kp];
117 Double_t pz = p * EMFSTK.wemf[kp];
118 Double_t tof = EMFSTK.agemf[kp];
119 Double_t polx = EMFSTK.upol[kp];
120 Double_t poly = EMFSTK.vpol[kp];
121 Double_t polz = EMFSTK.wpol[kp];
122 Double_t vx = EMFSTK.xemf[kp];
123 Double_t vy = EMFSTK.yemf[kp];
124 Double_t vz = EMFSTK.zemf[kp];
125 Double_t weight = EMFSTK.wtemf[kp];
131 //* case of no parent left (pair, photoelectric, annihilation):
132 //* all secondaries are true
133 if ((EVTFLG.lpairp == 1) || (EVTFLG.lphoel == 1) ||
134 (EVTFLG.lannfl == 1) || (EVTFLG.lannrs == 1)) {
136 if (EVTFLG.lpairp == 1) mech = kPPair;
137 else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
138 else mech = kPAnnihilation;
139 cppstack->PushTrack(done, parent, pdg,
140 px, py, pz, e, vx, vy, vz, tof,
141 polx, poly, polz, mech, ntr, weight, is);
142 if (debug) cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << " energy " << e-PAPROP.am[flukaid+6] << endl;
144 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
145 EMFSTK.iespak[kp][mkbmx2-2] = 0;
146 } // end of lpairp, lphoel, lannfl, lannrs
148 //* Compton: secondary is true only if charged (e+, e-)
149 else if ((EVTFLG.lcmptn == 1)) {
151 if (EMFSTK.ichemf[kp] != 0) {
153 cppstack->PushTrack(done, parent, pdg,
154 px, py, pz, e, vx, vy, vz, tof,
155 polx, poly, polz, mech, ntr, weight, is);
156 if (debug) cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
157 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
158 EMFSTK.iespak[kp][mkbmx2-2] = 0;
160 fluka->SetPint(px, py, pz, e);
164 //* Bremsstrahlung: true secondary only if charge = 0 (photon)
165 else if ((EVTFLG.lbrmsp == 1)) {
166 if (EMFSTK.ichemf[kp] == 0) {
168 cppstack->PushTrack(done, parent, pdg,
169 px, py, pz, e, vx, vy, vz, tof,
170 polx, poly, polz, mech, ntr, weight, is);
171 if (debug) cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
172 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
173 EMFSTK.iespak[kp][mkbmx2-2] = 0;
175 fluka->SetPint(px, py, pz, e);
180 //* Delta ray: If Bhabha, true secondary only if negative (electron)
181 else if ((EVTFLG.ldltry == 1)) {
183 if (EMFSTK.ichemf[kp] == -1) {
185 cppstack->PushTrack(done, parent, pdg,
186 px, py, pz, e, vx, vy, vz, tof,
187 polx, poly, polz, mech, ntr, weight, is);
188 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
189 EMFSTK.iespak[kp][mkbmx2-2] = 0;
190 if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
192 fluka->SetPint(px, py, pz, e);
196 //* Delta ray: Otherwise Moller: true secondary is the electron with
197 //* lower energy, which has been put higher in the stack
198 else if (kp == EMFSTK.npemf-1) {
200 cppstack->PushTrack(done, parent, pdg,
201 px, py, pz, e, vx, vy, vz, tof,
202 polx, poly, polz, mech, ntr, weight, is);
203 if (debug) cout << endl << " !!! stupre (Moller) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
204 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
205 EMFSTK.iespak[kp][mkbmx2-2] = 0;
207 fluka->SetPint(px, py, pz, e);
213 } // end of extern "C"