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