Use fluka routine usrdci to get ion properties.
[u/mrichter/AliRoot.git] / TFluka / stuprf.cxx
CommitLineData
f827183d 1#include <Riostream.h>
97b67dd4 2#include "TCallf77.h" //For the fortran calls
f827183d 3#ifndef WIN32
4# define stuprf stuprf_
97b67dd4 5# define usrdci usrdci_
f827183d 6#else
7# define stuprf STUPRF
97b67dd4 8# define usrdci USRDCI
f827183d 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
81f1d030 17#include "Fflkstk.h" //(FLKSTK) fluka common
f827183d 18#include "Ftrackr.h" //(TRACKR) fluka common
81f1d030 19#include "Fgenstk.h" //(GENSTK) fluka common
ca01d0af 20#include "Ffheavy.h" //(FHEAVY) fluka common
f827183d 21
18e0cabb 22
f827183d 23//Virtual MC
24#include "TFluka.h"
25#include "TVirtualMCStack.h"
26#include "TVirtualMCApplication.h"
27#include "TParticle.h"
28#include "TVector3.h"
29
30extern "C" {
97b67dd4 31
32 void type_of_call usrdci(int&, int&, int&, int&);
33
70541a80 34 void stuprf(Int_t& /*ij*/, Int_t& /*mreg*/,
18e0cabb 35 Double_t& xx, Double_t& yy, Double_t& zz,
36 Int_t& numsec, Int_t& npprmr)
f827183d 37{
38//*----------------------------------------------------------------------*
39//* *
40//* SeT User PRoperties for Fluka particles *
41//* *
42//*----------------------------------------------------------------------*
ca01d0af 43// Get the pointer to the VMC
44 TFluka* fluka = (TFluka*) gMC;
f827183d 45
81f1d030 46// FLKSTK.npflka = stack pointer
47// FLKSTK.louse = user flag
f827183d 48// TRACKR.llouse = user defined flag for the current particle
5834fcef 49
50 FLKSTK.louse[FLKSTK.npflka] = TRACKR.llouse;
f827183d 51
52// mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR
53// mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR
81f1d030 54// FLKSTK.sparek = spare real variables available for k.w.burn
55// FLKSTK.ispark = spare integer variables available for k.w.burn
f827183d 56// TRACKR.spausr = user defined spare variables for the current particle
57// TRACKR.ispusr = user defined spare flags for the current particle
58 Int_t ispr;
5834fcef 59 if (numsec <= npprmr) {
60 for (ispr = 0; ispr <= mkbmx1 - 1; ispr++) {
61 FLKSTK.sparek[FLKSTK.npflka][ispr] = TRACKR.spausr[ispr];
5834fcef 62 }
63 for (ispr = 0; ispr <= mkbmx2 - 1; ispr++) {
64 FLKSTK.ispark[FLKSTK.npflka][ispr] = TRACKR.ispusr[ispr];
5834fcef 65 }
66 }
67
18e0cabb 68 // save parent info
5834fcef 69 FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 3] = TRACKR.jtrack; // fluka particle id
18e0cabb 70 FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 4] = TRACKR.ispusr[mkbmx2 - 1]; // current track number
5834fcef 71 FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 5] = npprmr; // flag npprmr>0
18e0cabb 72
ca01d0af 73
2bc4c610 74 Int_t verbosityLevel = fluka->GetVerbosityLevel();
5834fcef 75 Bool_t debug = (verbosityLevel>=3)? kTRUE:kFALSE;
15b1e45e 76
fbf08100 77 fluka->SetTrackIsNew(kTRUE);
78// TVirtualMC* fluka = TFluka::GetMC();
f827183d 79// Get the stack produced from the generator
80 TVirtualMCStack* cppstack = fluka->GetStack();
81
82// EVTFLG.ntrcks = track number
83// Increment the track number and put it into the last flag
81ee7392 84// was numsec -1
85// clarify with Alberto
18e0cabb 86
87// npprmr > 0, the secondary being loaded is actually still the interacting
88// particle (it can happen in some biasing situations)
89
15b1e45e 90 if (numsec > npprmr) {
642f15cf 91// Now call the PushTrack(...)
70541a80 92 Int_t done = 0;
81ee7392 93
41737358 94 Int_t parent = TRACKR.ispusr[mkbmx2-1];
5834fcef 95 Int_t kpart = GENSTK.kpart[numsec-1];
ca01d0af 96 Int_t pdg;
97b67dd4 97 Int_t a = 1;
98 Int_t z = 1;
99 Int_t ism = 0;
ca01d0af 100 if (kpart < -6) {
97b67dd4 101// kpart = -kpart;
102// z = kpart / 100000;
103// a = kpart - z * 100000;
104// a /= 100;
105 usrdci(kpart, a, z, ism);
ca01d0af 106 pdg = fluka->GetIonPdg(z, a);
107 fluka->AddIon(a, z);
108 } else {
109 pdg = fluka->PDGFromId(kpart);
110 }
111
81f1d030 112 Double_t px = GENSTK.plr[numsec-1] * GENSTK.cxr[numsec-1];
113 Double_t py = GENSTK.plr[numsec-1] * GENSTK.cyr[numsec-1];
114 Double_t pz = GENSTK.plr[numsec-1] * GENSTK.czr[numsec-1];
ca01d0af 115 Double_t e;
116
97b67dd4 117 if (a == 1) {
ca01d0af 118 e = GENSTK.tki[numsec-1] + PAPROP.am[GENSTK.kpart[numsec-1]+6];
119 } else {
120 e = GENSTK.tki[numsec-1] + Double_t(a) * 0.93149 + 8.071e-3;
121 }
122
123
f827183d 124 Double_t vx = xx;
125 Double_t vy = yy;
126 Double_t vz = zz;
18e0cabb 127
81ee7392 128 Double_t tof = TRACKR.atrack;
81f1d030 129 Double_t polx = GENSTK.cxrpol[numsec-1];
130 Double_t poly = GENSTK.cyrpol[numsec-1];
131 Double_t polz = GENSTK.czrpol[numsec-1];
15b1e45e 132
f827183d 133 TMCProcess mech = kPHadronic;
81ee7392 134 if (EVTFLG.ldecay == 1) {
18e0cabb 135 mech = kPDecay;
136 if (debug) cout << endl << "Decay" << endl;
4aba9d66 137 FLKSTK.nlattc[FLKSTK.npflka] = TRACKR.lt1trk;
81ee7392 138 } else if (EVTFLG.ldltry == 1) {
18e0cabb 139 mech = kPDeltaRay;
140 if( fluka->GetIcode() == kKASHEA ) {
141 // For all interactions secondaries are put on GENSTK common (kp=1,np)
142 // but for KASHEA delta ray generation where only the secondary elec-
143 // tron is present and stacked on FLKSTK common for kp=lstack
144 pdg = fluka->PDGFromId( FLKSTK.iloflk[FLKSTK.npflka] );
145 px = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.txflk[FLKSTK.npflka];
146 py = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.tyflk[FLKSTK.npflka];
147 pz = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.tzflk[FLKSTK.npflka];
148 e = FLKSTK.tkeflk[FLKSTK.npflka] + PAPROP.am[FLKSTK.iloflk[FLKSTK.npflka]+6];
149 polx = FLKSTK.txpol[FLKSTK.npflka];
150 poly = FLKSTK.typol[FLKSTK.npflka];
151 polz = FLKSTK.tzpol[FLKSTK.npflka];
152 if (debug) cout << endl << "Delta Ray from KASHEA...." << " pdg from FLKSTK=" << pdg << endl;
15b1e45e 153 } else {
154 if (debug) cout << endl << "Delta Ray" << endl;
155 }
81ee7392 156 } else if (EVTFLG.lpairp == 1) {
18e0cabb 157 mech = kPPair;
158 if (debug) cout << endl << "Pair Production" << endl;
81ee7392 159 } else if (EVTFLG.lbrmsp == 1) {
18e0cabb 160 mech = kPBrem;
161 if (debug) cout << endl << "Bremsstrahlung" << endl;
81ee7392 162 }
81ee7392 163
81f1d030 164 Double_t weight = GENSTK.wei[numsec-1];
f827183d 165 Int_t is = 0;
18e0cabb 166 Int_t ntr;
81ee7392 167 //
168 // Save particle in VMC stack
642f15cf 169 cppstack->PushTrack(done, parent, pdg,
15b1e45e 170 px, py, pz, e,
171 vx, vy, vz, tof,
172 polx, poly, polz,
173 mech, ntr, weight, is);
ca01d0af 174 if (debug) {
175
15b1e45e 176 cout << endl << " !!! stuprf: ntr=" << ntr << " pdg " << pdg << " parent=" << parent
177 << " parent_pdg="<< fluka->PDGFromId(TRACKR.jtrack) << " numsec "
5834fcef 178 << numsec << " npprmr " << npprmr << " icode=" << fluka->GetIcode()
179 << "kin. energy [keV] = " << GENSTK.tki[numsec-1] * 1.e6 << endl << endl;
15b1e45e 180
ca01d0af 181 }
182
18e0cabb 183
81ee7392 184//
185// Save current track number
81f1d030 186 FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1] = ntr;
187 FLKSTK.ispark[FLKSTK.npflka][mkbmx2-2] = 0;
18e0cabb 188 } // end of if (numsec > npprmr)
189// else {
190// if(debug) {
191// cout << endl << " !!! stuprf: skipping pushtrack track=" << TRACKR.ispusr[mkbmx2-1]
192// << " pdg " << fluka->PDGFromId(TRACKR.jtrack) << " numsec=" << numsec<< " npprmr=" << npprmr
193// << " GENSTK pdg=" << fluka->PDGFromId(GENSTK.kpart[numsec-1]) << endl;
194// }
195// }
f827183d 196} // end of stuprf
97b67dd4 197
f827183d 198} // end of extern "C"
199