Parent id and current track id handling corrected.
[u/mrichter/AliRoot.git] / TFluka / stuprf.cxx
1 #include <Riostream.h>
2 #include "AliRun.h"
3 #include "AliStack.h"
4 #include "TFluka.h"
5 #ifndef WIN32
6 # define stuprf stuprf_
7 #else
8 # define stuprf STUPRF
9 #endif
10 //
11 // Fluka include
12 #include "Fdimpar.h"  //(DIMPAR) fluka include
13 // Fluka commons
14 #include "Fdblprc.h"  //(DBLPRC) fluka common
15 #include "Fevtflg.h"  //(EVTFLG) fluka common
16 #include "Fpaprop.h"  //(PAPROP) fluka common
17 #include "Fstack.h"   //(STACK)  fluka common
18 #include "Ftrackr.h"  //(TRACKR) fluka common
19 #include "Ffinuc.h"   //(FINUC)  fluka common
20
21 //Virtual MC
22 #include "TFluka.h"
23 #include "TVirtualMCStack.h"
24 #include "TVirtualMCApplication.h"
25 #include "TParticle.h"
26 #include "TVector3.h"
27
28 extern "C" {
29 void stuprf(Int_t& ij, Int_t& mreg,
30             Double_t& xx, Double_t& yy, Double_t& zz,
31             Int_t& numsec, Int_t& npprmr)
32 {
33 //*----------------------------------------------------------------------*
34 //*                                                                      *
35 //*  SeT User PRoperties for Fluka particles                             *
36 //*                                                                      *
37 //*----------------------------------------------------------------------*
38
39 // STACK.lstack  = stack pointer
40 // STACK.louse   = user flag
41 // TRACKR.llouse = user defined flag for the current particle
42   STACK.louse[STACK.lstack] = TRACKR.llouse;
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 // STACK.sparek  = spare real variables available for k.w.burn
47 // STACK.ispark  = spare integer variables available for k.w.burn
48 // TRACKR.spausr = user defined spare variables for the current particle
49 // TRACKR.ispusr = user defined spare flags for the current particle
50   Int_t ispr;
51   for (ispr = 0; ispr <= mkbmx1 - 1; ispr++) {
52     STACK.sparek[STACK.lstack][ispr] = TRACKR.spausr[ispr];
53   }  
54   for (ispr = 0; ispr <= mkbmx2 - 1; ispr++) {
55     STACK.ispark[STACK.lstack][ispr] = TRACKR.ispusr[ispr];
56   }  
57  
58 // Get the pointer to the VMC
59   TVirtualMC* fluka = TFluka::GetMC();
60 // Get the stack produced from the generator
61   TVirtualMCStack* cppstack = fluka->GetStack();
62   
63 // EVTFLG.ntrcks = track number
64 // Increment the track number and put it into the last flag
65 // was numsec -1
66 // clarify with Alberto
67   if (numsec > npprmr) {
68 // Now call the SetTrack(...)
69     Int_t done = 1;
70
71     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
72
73     Int_t pdg = fluka->PDGFromId(FINUC.kpart[numsec-1]);
74     Double_t px = FINUC.plr[numsec-1] * FINUC.cxr[numsec-1];
75     Double_t pz = FINUC.plr[numsec-1] * FINUC.cyr[numsec-1];
76     Double_t py = FINUC.plr[numsec-1] * FINUC.czr[numsec-1];
77     Double_t e  = FINUC.tki[numsec-1] + PAPROP.am[FINUC.kpart[numsec-1]+6];
78     Double_t vx = xx;
79     Double_t vy = yy;
80     Double_t vz = zz;
81     
82     Double_t tof  = TRACKR.atrack;
83     Double_t polx = FINUC.cxrpol[numsec-1];
84     Double_t poly = FINUC.cyrpol[numsec-1];
85     Double_t polz = FINUC.czrpol[numsec-1];
86     
87
88     TMCProcess mech = kPHadronic;
89     
90     if (EVTFLG.ldecay == 1) {
91         mech = kPDecay;
92         cout << endl << "Decay" << endl;
93         
94     } else if (EVTFLG.ldltry == 1) {
95         mech = kPDeltaRay;
96         cout << endl << "Delta Ray" << endl;
97         
98     } else if (EVTFLG.lpairp == 1) {
99         mech = kPPair;
100         cout << endl << "Pair Production" << endl;
101         
102     } else if (EVTFLG.lbrmsp == 1) {
103         mech = kPBrem;
104         cout << endl << "Bremsstrahlung" << endl;
105         
106     }
107     
108
109     Double_t weight = FINUC.wei[numsec-1];
110     Int_t is = 0;
111     Int_t ntr;  
112     // 
113     // Save particle in VMC stack
114     cppstack->SetTrack(done, parent, pdg,
115                        px, py, pz, e,
116                        vx, vy, vz, tof,
117                        polx, poly, polz,
118                        mech, ntr, weight, is);
119     cout << endl << " !!! stuprf: ntr=" << ntr << "pdg " << pdg << " parent=" << parent << "numsec " 
120          << numsec << "npprmr " << npprmr << endl;
121 //
122 //  Save current track number
123     STACK.ispark[STACK.lstack][mkbmx2-1] = ntr;
124   } // end of if (numsec-1 > npprmr)
125 } // end of stuprf
126 } // end of extern "C"
127