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