d876b9a6c9dd6024b0beaf056401cb5883e7859d
[u/mrichter/AliRoot.git] / TFluka / stuprf.cxx
1 #include <Riostream.h>
2 #include "AliRun.h"
3 #include "AliStack.h"
4 #ifndef WIN32
5 # define stuprf stuprf_
6 #else
7 # define stuprf STUPRF
8 #endif
9 //
10 // Fluka include
11 #include "Fdimpar.h"  //(DIMPAR) fluka include
12 // Fluka commons
13 #include "Fdblprc.h"  //(DBLPRC) fluka common
14 #include "Fevtflg.h"  //(EVTFLG) fluka common
15 #include "Fpaprop.h"  //(PAPROP) fluka common
16 #include "Fstack.h"   //(STACK)  fluka common
17 #include "Ftrackr.h"  //(TRACKR) fluka common
18 #include "Ffinuc.h"   //(FINUC)  fluka common
19
20 //Virtual MC
21 #include "TFluka.h"
22
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   TFluka* fluka =  (TFluka*) gMC;
60   Int_t verbosityLevel = fluka->GetVerbosityLevel();
61   Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
62   fluka->SetTrackIsNew(kTRUE);
63 //  TVirtualMC* fluka = TFluka::GetMC();
64 // Get the stack produced from the generator
65   TVirtualMCStack* cppstack = fluka->GetStack();
66   
67 // EVTFLG.ntrcks = track number
68 // Increment the track number and put it into the last flag
69 // was numsec -1
70 // clarify with Alberto
71   if (numsec > npprmr) {
72 // Now call the PushTrack(...)
73     Int_t done = 0;
74
75     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
76     Int_t kpart  = FINUC.kpart[numsec-1];
77     if (kpart < -6) return;
78
79     Int_t pdg = fluka->PDGFromId(kpart);
80
81     
82     Double_t px = FINUC.plr[numsec-1] * FINUC.cxr[numsec-1];
83     Double_t pz = FINUC.plr[numsec-1] * FINUC.cyr[numsec-1];
84     Double_t py = FINUC.plr[numsec-1] * FINUC.czr[numsec-1];
85     Double_t e  = FINUC.tki[numsec-1] + PAPROP.am[FINUC.kpart[numsec-1]+6];
86
87     Double_t vx = xx;
88     Double_t vy = yy;
89     Double_t vz = zz;
90     
91     Double_t tof  = TRACKR.atrack;
92     Double_t polx = FINUC.cxrpol[numsec-1];
93     Double_t poly = FINUC.cyrpol[numsec-1];
94     Double_t polz = FINUC.czrpol[numsec-1];
95     
96
97     TMCProcess mech = kPHadronic;
98     
99     if (EVTFLG.ldecay == 1) {
100         mech = kPDecay;
101         if (debug) cout << endl << "Decay" << endl;
102         
103     } else if (EVTFLG.ldltry == 1) {
104         mech = kPDeltaRay;
105         if (debug) cout << endl << "Delta Ray" << endl;
106         
107     } else if (EVTFLG.lpairp == 1) {
108         mech = kPPair;
109         if (debug) cout << endl << "Pair Production" << endl;
110         
111     } else if (EVTFLG.lbrmsp == 1) {
112         mech = kPBrem;
113         if (debug) cout << endl << "Bremsstrahlung" << endl;
114         
115     }
116     
117
118     Double_t weight = FINUC.wei[numsec-1];
119     Int_t is = 0;
120     Int_t ntr;  
121     // 
122     // Save particle in VMC stack
123     cppstack->PushTrack(done, parent, pdg,
124                        px, py, pz, e,
125                        vx, vy, vz, tof,
126                        polx, poly, polz,
127                        mech, ntr, weight, is);
128     if (debug) cout << endl << " !!! stuprf: ntr=" << ntr << "pdg " << pdg << " parent=" << parent << "numsec " 
129          << numsec << "npprmr " << npprmr << endl;
130 //
131 //  Save current track number
132     STACK.ispark[STACK.lstack][mkbmx2-1] = ntr;
133     STACK.ispark[STACK.lstack][mkbmx2-2] = 0;
134   } // end of if (numsec-1 > npprmr)
135 } // end of stuprf
136 } // end of extern "C"
137