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