Moving to the new VMC naming convention
[u/mrichter/AliRoot.git] / TFluka / stupre.cxx
1 #include <Riostream.h>
2 #include "AliRun.h"
3 #include "TFluka.h"
4 #ifndef WIN32
5 # define stupre stupre_
6 #else
7 # define stupre STUPRE
8 #endif
9 //
10 // Fluka include
11 #include "Fdimpar.h"  //(DIMPAR) fluka include
12 // Fluka commons
13 #include "Fdblprc.h"  //(DBLPRC) fluka common
14 #include "Femfstk.h"  //(EMFSTK) fluka common
15 #include "Fevtflg.h"  //(EVTFLG) fluka common
16 #include "Fpaprop.h"  //(PAPROP) fluka common
17 #include "Ftrackr.h"  //(TRACKR) fluka common
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
26 extern "C" {
27 void stupre()
28 {
29 //*----------------------------------------------------------------------*
30 //*                                                                      *
31 //*  SeT User PRoperties for Emf particles                               *
32 //*                                                                      *
33 //*----------------------------------------------------------------------*
34
35   Int_t lbhabh = 0;
36   if (EVTFLG.ldltry == 1) {
37     if (EMFSTK.ichemf[EMFSTK.npemf-1] * EMFSTK.ichemf[EMFSTK.npemf-2] < 0) lbhabh = 1;
38   }
39
40 // mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR
41 // mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR
42 // EMFSTK.espark  = spare real variables available for 
43 // EMFSTK.iespak  = spare integer variables available for
44 // TRACKR.spausr = user defined spare variables for the current particle
45 // TRACKR.ispusr = user defined spare flags for the current particle
46 // EMFSTK.louemf = user flag
47 // TRACKR.llouse = user defined flag for the current particle
48
49   Int_t npnw, ispr;
50   for (npnw=EMFSTK.npstrt-1; npnw<=EMFSTK.npemf-1; npnw++) {
51
52     for (ispr=0; ispr<=mkbmx1-1; ispr++) 
53       EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr];
54
55     for (ispr=0; ispr<=mkbmx2-1; ispr++) 
56       EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr];
57
58     EMFSTK.louemf[npnw] = TRACKR.llouse;
59   }
60
61 // Get the pointer to the VMC
62   TVirtualMC* fluka = TFluka::GetMC();
63 // Get the stack produced from the generator
64   TVirtualMCStack* cppstack = fluka->GetStack();
65   
66 // EVTFLG.ntrcks = track number
67 // Increment the track number and put it into the last flag
68
69   Int_t kp;
70   for (kp = EMFSTK.npstrt - 1; kp <= EMFSTK.npemf - 1; kp++) {
71
72 //* save the parent track number and reset it at each loop
73     Int_t done = 1;
74
75     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
76     
77     Int_t flukaid = 0;
78
79     if (EMFSTK.ichemf[kp] == -1) flukaid = 3;
80     else if (EMFSTK.ichemf[kp] == 0)  flukaid = 7;
81     else if (EMFSTK.ichemf[kp] == 0)  flukaid = 4;
82     Int_t pdg       = fluka->PDGFromId(flukaid);
83     Double_t e      = EMFSTK.etemf[kp] * emvgev;
84     Double_t p      = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
85     Double_t px     = p * EMFSTK.u[kp];
86     Double_t pz     = p * EMFSTK.v[kp];
87     Double_t py     = p * EMFSTK.w[kp];
88     Double_t tof    = EMFSTK.agemf[kp];
89     Double_t polx   = EMFSTK.upol[kp];
90     Double_t poly   = EMFSTK.vpol[kp];
91     Double_t polz   = EMFSTK.wpol[kp];
92     Double_t vx     = EMFSTK.x[kp];
93     Double_t vy     = EMFSTK.y[kp];
94     Double_t vz     = EMFSTK.z[kp];
95     Double_t weight = EMFSTK.wtemf[kp];
96
97     Int_t ntr;
98     TMCProcess mech;
99     Int_t is = 0;
100
101 //* case of no parent left (pair, photoelectric, annihilation):
102 //* all secondaries are true
103     if ((EVTFLG.lpairp == 1) || (EVTFLG.lphoel == 1) ||
104         (EVTFLG.lannfl == 1) || (EVTFLG.lannrs == 1)) {
105         
106         if (EVTFLG.lpairp == 1) mech = kPPair;
107         else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
108         else mech = kPAnnihilation;
109         cppstack->PushTrack(done, parent, pdg,
110                            px, py, pz, e, vx, vy, vz, tof,
111                            polx, poly, polz, mech, ntr, weight, is);
112         cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
113
114         EMFSTK.iespak[kp][mkbmx2-1] = ntr;
115     } // end of lpairp, lphoel, lannfl, lannrs
116     
117 //* Compton: secondary is true only if charged (e+, e-)
118     else if ((EVTFLG.lcmptn == 1)) {
119
120         if (EMFSTK.ichemf[kp] != 0) {
121             mech = kPCompton;
122             cppstack->PushTrack(done, parent, pdg,
123                                px, py, pz, e, vx, vy, vz, tof,
124                                polx, poly, polz, mech, ntr, weight, is);
125             cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
126             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
127         }
128     } // end of lcmptn
129     
130 //* Bremsstrahlung: true secondary only if charge = 0 (photon)
131     else if ((EVTFLG.lbrmsp == 1)) {
132         if (EMFSTK.ichemf[kp] == 0) {
133             mech = kPBrem;
134             cppstack->PushTrack(done, parent, pdg,
135                                px, py, pz, e, vx, vy, vz, tof,
136                                polx, poly, polz, mech, ntr, weight, is);
137             cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
138             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
139         }
140     } // end of lbrmsp
141     
142 //* Delta ray: If Bhabha, true secondary only if negative (electron)
143     else if ((EVTFLG.ldltry == 1)) {
144         if (lbhabh == 1) {
145             if (EMFSTK.ichemf[kp] == -1) {
146                 mech = kPDeltaRay;
147                 cppstack->PushTrack(done, parent, pdg,
148                                    px, py, pz, e, vx, vy, vz, tof,
149                                    polx, poly, polz, mech, ntr, weight, is);
150                 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
151             cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
152             } // end of Bhabha
153         } // lbhabh == 1
154         
155 //* Delta ray: Otherwise Moller: true secondary is the electron with
156 //*            lower energy, which has been put higher in the stack
157         else if (kp == EMFSTK.npemf-1) {
158             mech = kPDeltaRay;
159             cppstack->PushTrack(done, parent, pdg,
160                                px, py, pz, e, vx, vy, vz, tof,
161                                polx, poly, polz, mech, ntr, weight, is);
162             cout << endl << " !!! stupre (Moller) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
163             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
164         } // end of Delta ray
165     } // end of ldltry
166     
167   } // end of loop
168   
169 // !!! TO BE CONFIRMED !!!
170 } // end of stupre
171 } // end of extern "C"
172