Update to FLUKA 2006.3b.
[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   fluka->SetTrackIsNew(kTRUE);
69
70 // Get the stack produced from the generator
71   TVirtualMCStack* cppstack = fluka->GetStack();
72   
73 // EVTFLG.ntrcks = track number
74 // Increment the track number and put it into the last flag
75
76   Int_t kp;
77   
78   for (kp = EMFSTK.npstrt - 1; kp <= EMFSTK.npemf - 1; kp++) {
79     
80 // Check transport cut first
81     Int_t    ireg   = EMFSTK.iremf[kp];
82     Double_t cut    = (TMath::Abs(EMFSTK.ichemf[kp]) == 1) ? EMFRGN.elethr[ireg-1] :  EMFRGN.phothr[ireg-1];
83     Double_t e      = EMFSTK.etemf[kp];
84
85     if ((e < cut) 
86         && (
87             (EMFSTK.ichemf[kp] ==  0) ||
88             (EMFSTK.ichemf[kp] == -1) ||
89             (EMFSTK.ichemf[kp] ==  1 &&  EMFRGN.phothr[ireg-1] > emassmev)
90             )
91         )
92     {
93         EMFSTK.iespak[kp][mkbmx2-1] = TRACKR.ispusr[mkbmx2-1];
94         EMFSTK.iespak[kp][mkbmx2-2] =  0;
95         EMFSTK.iespak[kp][mkbmx2-3] = TRACKR.jtrack;
96         EMFSTK.irlatt[kp] = TRACKR.lt1trk;
97         continue;
98     }
99
100 // Save the parent track number and reset it at each loop
101     Int_t done = 0;
102     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
103     Int_t flukaid = 0;
104
105 // Identify particle type
106     if      (EMFSTK.ichemf[kp] == -1)  flukaid = kFLUKAelectron;
107     else if (EMFSTK.ichemf[kp] ==  0)  flukaid = kFLUKAphoton;
108     else if (EMFSTK.ichemf[kp] ==  1)  flukaid = kFLUKApositron;
109
110
111     e *= emvgev;
112     Int_t    pdg    = fluka->PDGFromId(flukaid);
113     Double_t p      = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
114     Double_t px     = p * EMFSTK.uemf[kp];
115     Double_t py     = p * EMFSTK.vemf[kp];
116     Double_t pz     = p * EMFSTK.wemf[kp];
117     Double_t tof    = EMFSTK.agemf[kp];
118     Double_t polx   = EMFSTK.upol[kp];
119     Double_t poly   = EMFSTK.vpol[kp];
120     Double_t polz   = EMFSTK.wpol[kp];
121     Double_t vx     = EMFSTK.xemf[kp];
122     Double_t vy     = EMFSTK.yemf[kp];
123     Double_t vz     = EMFSTK.zemf[kp];
124     Double_t weight = EMFSTK.wtemf[kp];
125
126     Int_t ntr;
127     TMCProcess mech;
128     Int_t is = 0;
129
130 //* case of no parent left (pair, photoelectric, annihilation):
131 //* all secondaries are true
132     if ((EVTFLG.lpairp == 1) || (EVTFLG.lphoel == 1) ||
133         (EVTFLG.lannfl == 1) || (EVTFLG.lannrs == 1)) {
134
135         if (EVTFLG.lpairp == 1) mech = kPPair;
136         else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
137         else mech = kPAnnihilation;
138         cppstack->PushTrack(done, parent, pdg,
139                            px, py, pz, e, vx, vy, vz, tof,
140                            polx, poly, polz, mech, ntr, weight, is);
141         if (debug) cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << " energy " << e-PAPROP.am[flukaid+6] << endl;
142
143         EMFSTK.iespak[kp][mkbmx2-1] = ntr;
144         EMFSTK.iespak[kp][mkbmx2-2] = 0;
145     } // end of lpairp, lphoel, lannfl, lannrs
146
147 //* Compton: secondary is true only if charged (e+, e-)
148     else if ((EVTFLG.lcmptn == 1)) {
149
150         if (EMFSTK.ichemf[kp] != 0) {
151             mech = kPCompton;
152             cppstack->PushTrack(done, parent, pdg,
153                                px, py, pz, e, vx, vy, vz, tof,
154                                polx, poly, polz, mech, ntr, weight, is);
155             if (debug) cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
156             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
157             EMFSTK.iespak[kp][mkbmx2-2] = 0;
158         }
159     } // end of lcmptn
160
161 //* Bremsstrahlung: true secondary only if charge = 0 (photon)
162     else if ((EVTFLG.lbrmsp == 1)) {
163         if (EMFSTK.ichemf[kp] == 0) {
164             mech = kPBrem;
165             cppstack->PushTrack(done, parent, pdg,
166                                px, py, pz, e, vx, vy, vz, tof,
167                                polx, poly, polz, mech, ntr, weight, is);
168             if (debug) cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
169             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
170             EMFSTK.iespak[kp][mkbmx2-2] = 0;
171         }
172     } // end of lbrmsp
173
174 //* Delta ray: If Bhabha, true secondary only if negative (electron)
175     else if ((EVTFLG.ldltry == 1)) {
176         if (lbhabh == 1) {
177             if (EMFSTK.ichemf[kp] == -1) {
178                 mech = kPDeltaRay;
179                 cppstack->PushTrack(done, parent, pdg,
180                                    px, py, pz, e, vx, vy, vz, tof,
181                                    polx, poly, polz, mech, ntr, weight, is);
182                 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
183                 EMFSTK.iespak[kp][mkbmx2-2] = 0;
184            if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
185             } // end of Bhabha
186         } // lbhabh == 1
187
188 //* Delta ray: Otherwise Moller: true secondary is the electron with
189 //*            lower energy, which has been put higher in the stack
190         else if (kp == EMFSTK.npemf-1) {
191             mech = kPDeltaRay;
192             cppstack->PushTrack(done, parent, pdg,
193                                px, py, pz, e, vx, vy, vz, tof,
194                                polx, poly, polz, mech, ntr, weight, is);
195             if (debug) cout << endl << " !!! stupre (Moller) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
196             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
197             EMFSTK.iespak[kp][mkbmx2-2] = 0;
198         } // end of Delta ray
199     } // end of ldltry
200
201   } // end of loop
202 } // end of stupre
203 } // end of extern "C"
204