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