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