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