Correct energy of interacting particle after BREMS, DELTA, COMPTON
[u/mrichter/AliRoot.git] / TFluka / stupre.cxx
1 #include <Riostream.h>
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
10
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
17 #include "Femfrgn.h"  //(EFMRGN) fluka common
18
19 //Virtual MC
20 #include "TFluka.h"
21 #include "TFlukaCodes.h"
22 #include "TVirtualMCStack.h"
23 #include "TVirtualMCApplication.h"
24 #include "TParticle.h"
25 #include "TVector3.h"
26
27 extern "C" {
28 void stupre()
29 {
30 //*----------------------------------------------------------------------*
31 //*                                                                      *
32 //*  SeT User PRoperties for Emf particles                               *
33 //*                                                                      *
34 //*----------------------------------------------------------------------*
35 // Get the pointer to the VMC
36   TFluka* fluka =  (TFluka*) gMC;
37   static Double_t emassmev = PAPROP.am[kFLUKAelectron + 6] * 1000.;
38     
39   Int_t lbhabh = 0;
40   if (EVTFLG.ldltry == 1) {
41     if (EMFSTK.ichemf[EMFSTK.npemf-1] * EMFSTK.ichemf[EMFSTK.npemf-2] < 0) lbhabh = 1;
42   }
43
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
52
53   Int_t npnw, ispr;
54   for (npnw = EMFSTK.npstrt-1; npnw <= EMFSTK.npemf-1; npnw++) {
55
56       for (ispr = 0; ispr <= mkbmx1-1; ispr++) 
57           EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr];
58
59       for (ispr = 0; ispr <= mkbmx2-1; ispr++) 
60           EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr];
61
62       EMFSTK.louemf[npnw] = TRACKR.llouse;
63   }
64
65
66   Int_t verbosityLevel = fluka->GetVerbosityLevel();
67   Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
68   debug = 100;
69   
70   fluka->SetTrackIsNew(kTRUE);
71
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;
79   
80   for (kp = EMFSTK.npstrt - 1; kp <= EMFSTK.npemf - 1; kp++) {
81     
82 // Check transport cut first
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];
85     Double_t e      = EMFSTK.etemf[kp];
86
87     if ((e < cut) 
88         && (
89             (EMFSTK.ichemf[kp] ==  0) ||
90             (EMFSTK.ichemf[kp] == -1) ||
91             (EMFSTK.ichemf[kp] ==  1 &&  EMFRGN.phothr[ireg-1] > emassmev)
92             )
93         )
94     {
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;
100     }
101
102 // Save the parent track number and reset it at each loop
103     Int_t done = 0;
104     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
105     Int_t flukaid = 0;
106
107 // Identify particle type
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;
111
112
113     e *= emvgev;
114     Int_t    pdg    = fluka->PDGFromId(flukaid);
115     Double_t p      = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
116     Double_t px     = p * EMFSTK.uemf[kp];
117     Double_t py     = p * EMFSTK.vemf[kp];
118     Double_t pz     = p * EMFSTK.wemf[kp];
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];
123     Double_t vx     = EMFSTK.xemf[kp];
124     Double_t vy     = EMFSTK.yemf[kp];
125     Double_t vz     = EMFSTK.zemf[kp];
126     Double_t weight = EMFSTK.wtemf[kp];
127
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)) {
136
137         if (EVTFLG.lpairp == 1) mech = kPPair;
138         else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
139         else mech = kPAnnihilation;
140         cppstack->PushTrack(done, parent, pdg,
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;
144
145         EMFSTK.iespak[kp][mkbmx2-1] = ntr;
146         EMFSTK.iespak[kp][mkbmx2-2] = 0;
147     } // end of lpairp, lphoel, lannfl, lannrs
148
149 //* Compton: secondary is true only if charged (e+, e-)
150     else if ((EVTFLG.lcmptn == 1)) {
151
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;
160         } else {
161             fluka->SetPint(px, py, pz, e);
162         }
163     } // end of lcmptn
164
165 //* Bremsstrahlung: true secondary only if charge = 0 (photon)
166     else if ((EVTFLG.lbrmsp == 1)) {
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;
175         } else {
176             fluka->SetPint(px, py, pz, e);
177         }
178         
179     } // end of lbrmsp
180
181 //* Delta ray: If Bhabha, true secondary only if negative (electron)
182     else if ((EVTFLG.ldltry == 1)) {
183         if (lbhabh == 1) {
184             if (EMFSTK.ichemf[kp] == -1) {
185                 mech = kPDeltaRay;
186                 cppstack->PushTrack(done, parent, pdg,
187                                     px, py, pz, e, vx, vy, vz, tof,
188                                     polx, poly, polz, mech, ntr, weight, is);
189                 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
190                 EMFSTK.iespak[kp][mkbmx2-2] = 0;
191                 if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
192             } else {
193                 fluka->SetPint(px, py, pz, e);
194             }
195         } // lbhabh == 1
196
197 //* Delta ray: Otherwise Moller: true secondary is the electron with
198 //*            lower energy, which has been put higher in the stack
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;
207         } else {
208             fluka->SetPint(px, py, pz, e);
209         }
210     } // end of ldltry
211
212   } // end of loop
213 } // end of stupre
214 } // end of extern "C"
215