Added preprocessor conditionals to support ROOT > 5.11.2.
[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 // Ckeck 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         continue;
97     }
98
99 // Save the parent track number and reset it at each loop
100     Int_t done = 0;
101     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
102     Int_t flukaid = 0;
103     
104 // Identify particle type
105     if      (EMFSTK.ichemf[kp] == -1)  flukaid = kFLUKAelectron;
106     else if (EMFSTK.ichemf[kp] ==  0)  flukaid = kFLUKAphoton;
107     else if (EMFSTK.ichemf[kp] ==  1)  flukaid = kFLUKApositron;
108     
109
110     e *= emvgev;
111     Int_t    pdg    = fluka->PDGFromId(flukaid);
112     Double_t p      = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
113     Double_t px     = p * EMFSTK.uemf[kp];
114     Double_t py     = p * EMFSTK.vemf[kp];
115     Double_t pz     = p * EMFSTK.wemf[kp];
116     Double_t tof    = EMFSTK.agemf[kp];
117     Double_t polx   = EMFSTK.upol[kp];
118     Double_t poly   = EMFSTK.vpol[kp];
119     Double_t polz   = EMFSTK.wpol[kp];
120     Double_t vx     = EMFSTK.xemf[kp];
121     Double_t vy     = EMFSTK.yemf[kp];
122     Double_t vz     = EMFSTK.zemf[kp];
123     Double_t weight = EMFSTK.wtemf[kp];
124
125     Int_t ntr;
126     TMCProcess mech;
127     Int_t is = 0;
128
129 //* case of no parent left (pair, photoelectric, annihilation):
130 //* all secondaries are true
131     if ((EVTFLG.lpairp == 1) || (EVTFLG.lphoel == 1) ||
132         (EVTFLG.lannfl == 1) || (EVTFLG.lannrs == 1)) {
133         
134         if (EVTFLG.lpairp == 1) mech = kPPair;
135         else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
136         else mech = kPAnnihilation;
137         cppstack->PushTrack(done, parent, pdg,
138                            px, py, pz, e, vx, vy, vz, tof,
139                            polx, poly, polz, mech, ntr, weight, is);
140         if (debug) cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << "energy " << e-PAPROP.am[flukaid+6] << endl;
141
142         EMFSTK.iespak[kp][mkbmx2-1] = ntr;
143         EMFSTK.iespak[kp][mkbmx2-2] = 0;
144     } // end of lpairp, lphoel, lannfl, lannrs
145     
146 //* Compton: secondary is true only if charged (e+, e-)
147     else if ((EVTFLG.lcmptn == 1)) {
148
149         if (EMFSTK.ichemf[kp] != 0) {
150             mech = kPCompton;
151             cppstack->PushTrack(done, parent, pdg,
152                                px, py, pz, e, vx, vy, vz, tof,
153                                polx, poly, polz, mech, ntr, weight, is);
154             if (debug) cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
155             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
156             EMFSTK.iespak[kp][mkbmx2-2] = 0;
157         }
158     } // end of lcmptn
159     
160 //* Bremsstrahlung: true secondary only if charge = 0 (photon)
161     else if ((EVTFLG.lbrmsp == 1)) {
162         if (EMFSTK.ichemf[kp] == 0) {
163             mech = kPBrem;
164             cppstack->PushTrack(done, parent, pdg,
165                                px, py, pz, e, vx, vy, vz, tof,
166                                polx, poly, polz, mech, ntr, weight, is);
167             if (debug) cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
168             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
169             EMFSTK.iespak[kp][mkbmx2-2] = 0;
170         }
171     } // end of lbrmsp
172     
173 //* Delta ray: If Bhabha, true secondary only if negative (electron)
174     else if ((EVTFLG.ldltry == 1)) {
175         if (lbhabh == 1) {
176             if (EMFSTK.ichemf[kp] == -1) {
177                 mech = kPDeltaRay;
178                 cppstack->PushTrack(done, parent, pdg,
179                                    px, py, pz, e, vx, vy, vz, tof,
180                                    polx, poly, polz, mech, ntr, weight, is);
181                 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
182                 EMFSTK.iespak[kp][mkbmx2-2] = 0;
183            if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
184             } // end of Bhabha
185         } // lbhabh == 1
186         
187 //* Delta ray: Otherwise Moller: true secondary is the electron with
188 //*            lower energy, which has been put higher in the stack
189         else if (kp == EMFSTK.npemf-1) {
190             mech = kPDeltaRay;
191             cppstack->PushTrack(done, parent, pdg,
192                                px, py, pz, e, vx, vy, vz, tof,
193                                polx, poly, polz, mech, ntr, weight, is);
194             if (debug) cout << endl << " !!! stupre (Moller) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
195             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
196             EMFSTK.iespak[kp][mkbmx2-2] = 0;
197         } // end of Delta ray
198     } // end of ldltry
199     
200   } // end of loop
201 } // end of stupre
202 } // end of extern "C"
203