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