Parent id and current track id handling corrected.
[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.iq[EMFSTK.np-1] * EMFSTK.iq[EMFSTK.np-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.np-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; npnw <= EMFSTK.np-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     if (EMFSTK.iq[kp] == -1) flukaid = 3;
79     else if (EMFSTK.iq[kp] == 0)  flukaid = 7;
80     else if (EMFSTK.iq[kp] == 0)  flukaid = 4;
81     Int_t pdg = fluka->PDGFromId(flukaid);
82     Double_t e  = EMFSTK.e[kp]*emvgev;
83     Double_t p = sqrt(e*e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
84     Double_t px = p * EMFSTK.u[kp];
85     Double_t pz = p * EMFSTK.v[kp];
86     Double_t py = p * EMFSTK.w[kp];
87     Double_t tof = EMFSTK.agemf[kp];
88     Double_t polx = EMFSTK.upol[kp];
89     Double_t poly = EMFSTK.vpol[kp];
90     Double_t polz = EMFSTK.wpol[kp];
91     Double_t vx = EMFSTK.x[kp];
92     Double_t vy = EMFSTK.y[kp];
93     Double_t vz = EMFSTK.z[kp];
94     Double_t weight = EMFSTK.wt[kp];
95     Int_t ntr;
96     TMCProcess mech;
97     Int_t is = 0;
98
99 //* case of no parent left (pair, photoelectric, annihilation):
100 //* all secondaries are true
101     if ((EVTFLG.lpairp == 1) || (EVTFLG.lphoel == 1) ||
102         (EVTFLG.lannfl == 1) || (EVTFLG.lannrs == 1)) {
103         
104         if (EVTFLG.lpairp == 1) mech = kPPair;
105         else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
106         else mech = kPAnnihilation;
107         cppstack->SetTrack(done, parent, pdg,
108                            px, py, pz, e, vx, vy, vz, tof,
109                            polx, poly, polz, mech, ntr, weight, is);
110         
111         cout << endl << " !!! stupre: ntr=" << ntr  << " parent=" << parent << endl;
112         EMFSTK.iespak[kp][mkbmx2-1] = ntr;
113     } // end of lpairp, lphoel, lannfl, lannrs
114     
115 //* Compton: secondary is true only if charged (e+, e-)
116     else if ((EVTFLG.lcmptn == 1)) {
117         if (EMFSTK.iq[kp] != 0) {
118             mech = kPCompton;
119             cppstack->SetTrack(done, parent, pdg,
120                                px, py, pz, e, vx, vy, vz, tof,
121                                polx, poly, polz, mech, ntr, weight, is);
122             cout << endl << " !!! stupre: ntr=" << ntr  << " parent=" << parent << endl;
123             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
124         }
125     } // end of lcmptn
126     
127 //* Bremsstrahlung: true secondary only if charge = 0 (photon)
128     else if ((EVTFLG.lbrmsp == 1)) {
129         if (EMFSTK.iq[kp] == 0) {
130             mech = kPBrem;
131             cppstack->SetTrack(done, parent, pdg,
132                                px, py, pz, e, vx, vy, vz, tof,
133                                polx, poly, polz, mech, ntr, weight, is);
134             cout << endl << " !!! stupre: ntr=" << ntr  << " parent=" << parent << endl;
135             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
136         }
137     } // end of lbrmsp
138     
139 //* Delta ray: If Bhabha, true secondary only if negative (electron)
140     else if ((EVTFLG.ldltry == 1)) {
141         if (lbhabh == 1) {
142             if (EMFSTK.iq[kp] == 0) {
143                 mech = kPDeltaRay;
144                 cppstack->SetTrack(done, parent, pdg,
145                                    px, py, pz, e, vx, vy, vz, tof,
146                                    polx, poly, polz, mech, ntr, weight, is);
147                 EMFSTK.iespak[kp][mkbmx2-1] = ntr;
148             } // end of Bhabha
149         }
150         
151 //* Delta ray: Otherwise Moller: true secondary is the electron with
152 //*            lower energy, which has been put higher in the stack
153         else if (kp == EMFSTK.np) {
154             mech = kPDeltaRay;
155             cppstack->SetTrack(done, parent, pdg,
156                                px, py, pz, e, vx, vy, vz, tof,
157                                polx, poly, polz, mech, ntr, weight, is);
158             cout << endl << " !!! stupre: ntr=" << ntr << " parent=" << parent << endl;
159             EMFSTK.iespak[kp][mkbmx2-1] = ntr;
160         } // end of Delta ray
161     } // end of ldltry
162     
163   } // end of loop
164   
165 // !!! TO BE CONFIRMED !!!
166   EVTFLG.ntrcks = EMFSTK.iespak[EMFSTK.np-1][mkbmx2-1];
167   
168 } // end of stupre
169 } // end of extern "C"
170