Debug print removed.
[u/mrichter/AliRoot.git] / TFluka / source.cxx
CommitLineData
b9d0a01d 1// Fortran
2#include "TCallf77.h"
3
4// Fluka commons
5#include "Fdblprc.h" //(DBLPRC) fluka common
6#include "Fdimpar.h" //(DIMPAR) fluka parameters
81f1d030 7#include "Fsourcm.h" //(EPISOR) fluka common
7c33fcd7 8#include "Fflkstk.h" //(FLKSTK) fluka common
9#include "Fsumcou.h" //(SUMCOU) fluka common
b9d0a01d 10#include "Fpaprop.h" //(PAPROP) fluka common
11#include "Fltclcm.h" //(LTCLCM) fluka common
cfd35035 12#include "Fopphst.h" //(OPPHST) fluka common
7c33fcd7 13#include "Fioiocm.h" //(IOIOCM) fluka common
14#include "Fbeamcm.h" //(BEAMCM) fluka common
b9d0a01d 15//Virutal MC
16#include "TFluka.h"
7c33fcd7 17#include "TFlukaIon.h"
a7bb59a2 18
b9d0a01d 19#include "TVirtualMCStack.h"
fbf08100 20//#include "TVirtualMCApplication.h"
fbf08100 21
b9d0a01d 22#include "TParticle.h"
23#include "TVector3.h"
24
25//Other
eae0fe66 26#include <Riostream.h>
b9d0a01d 27
28#ifndef WIN32
29# define source source_
30# define geocrs geocrs_
31# define georeg georeg_
32# define geohsm geohsm_
33# define soevsv soevsv_
7c33fcd7 34# define dcdion dcdion_
35# define setion setion_
5fc9f5d5 36# define stisbm stisbm_
b9d0a01d 37#else
38# define source SOURCE
39# define geocrs GEOCRS
40# define georeg GEOREG
41# define geohsm GEOHSM
42# define soevsv SOEVSV
7c33fcd7 43# define dcdion DCDION
44# define setion SETION
5fc9f5d5 45# define stisbm STISBM
b9d0a01d 46#endif
b9d0a01d 47extern "C" {
48 //
49 // Prototypes for FLUKA functions
50 //
51 void type_of_call geocrs(Double_t &, Double_t &, Double_t &);
52 void type_of_call georeg(Double_t &, Double_t &, Double_t &,
4aba9d66 53 Int_t &, Int_t &);
b9d0a01d 54 void type_of_call geohsm(Int_t &, Int_t &, Int_t &, Int_t &);
55 void type_of_call soevsv();
7c33fcd7 56 void type_of_call dcdion(Int_t &);
57 void type_of_call setion(Int_t &);
5fc9f5d5 58 void type_of_call stisbm(Int_t &, Int_t &, Int_t &);
7c33fcd7 59
b9d0a01d 60 /*
61 *----------------------------------------------------------------------*
62 * *
63 * Created on 07 january 1990 by Alfredo Ferrari & Paola Sala *
64 * Infn - Milan *
65 * *
66 * Last change on 21-jun-98 by Alfredo Ferrari *
67 * *
68 * C++ version on 27-sep-02 by Isidro Gonzalez *
69 * *
70 * This is just an example of a possible user written source routine. *
71 * note that the beam card still has some meaning - in the scoring the *
72 * maximum momentum used in deciding the binning is taken from the *
73 * beam momentum. Other beam card parameters are obsolete. *
74 * *
75 *----------------------------------------------------------------------*/
76
77 void source(Int_t& nomore) {
b496f27c 78// Get the pointer to TFluka
2bc4c610 79 TFluka* fluka = (TFluka*)gMC;
b496f27c 80
2bc4c610 81 Int_t verbosityLevel = fluka->GetVerbosityLevel();
82 Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
83 if (debug) {
ce60a136 84 cout << "==> source(" << nomore << ")" << endl;
81f1d030 85 cout << "\t* SOURCM.lsouit = " << (SOURCM.lsouit?'T':'F') << endl;
2bc4c610 86 }
b9d0a01d 87
2bc4c610 88 static Bool_t lfirst = true;
89 static Bool_t particleIsPrimary = true;
90 static Bool_t lastParticleWasPrimary = true;
b9d0a01d 91
ce60a136 92 nomore = 0;
cadd9e6b 93
b496f27c 94
95// Get the stack
b9d0a01d 96 TVirtualMCStack* cppstack = fluka->GetStack();
b496f27c 97
ce60a136 98 TParticle* particle;
b9d0a01d 99 Int_t itrack = -1;
ce60a136 100 Int_t nprim = cppstack->GetNprimary();
101// Get the next particle from the stack
102 particle = cppstack->PopNextTrack(itrack);
2bc4c610 103 fluka->SetTrackIsNew(kTRUE);
a32460fe 104 if (itrack == (nprim - 1)) lfirst = true;
ce60a136 105// Is this a secondary not handled by Fluka, i.e. a particle added by user action ?
106 lastParticleWasPrimary = particleIsPrimary;
107
108 if (itrack >= nprim) {
4aba9d66 109 particleIsPrimary = kFALSE;
ce60a136 110 } else {
4aba9d66 111 particleIsPrimary = kTRUE;
ce60a136 112 }
113
ce60a136 114 if (lfirst) {
4aba9d66 115 SOURCM.tkesum = zerzer;
116 lfirst = false;
117 SOURCM.lussrc = true;
ce60a136 118 } else {
119//
120// Post-track actions for primary track
121//
4aba9d66 122 if (particleIsPrimary) {
123 TVirtualMCApplication::Instance()->PostTrack();
124 TVirtualMCApplication::Instance()->FinishPrimary();
125 if ((itrack%10)==0)
126 cout << "=== TRACKING PRIMARY "<< itrack <<" ===" << endl;
127 //printf("=== TRACKING PRIMARY %d ===\n", itrack);
128 }
ce60a136 129 }
b9d0a01d 130
b496f27c 131 // Exit if itrack is negative (-1). Set lsouit to false to mark last track for this event
ce60a136 132
b9d0a01d 133 if (itrack<0) {
134 nomore = 1;
81f1d030 135 SOURCM.lsouit = false;
2bc4c610 136 if (debug) {
81f1d030 137 cout << "\t* SOURCM.lsouit = " << (SOURCM.lsouit?'T':'F') << endl;
2bc4c610 138 cout << "\t* No more particles. Exiting..." << endl;
139 cout << "<== source(" << nomore << ")" << endl;
140 }
b9d0a01d 141 return;
142 }
ce60a136 143
b496f27c 144 //
145 // Handle user event abortion
146 if (fluka->EventIsStopped()) {
4aba9d66 147 printf("Event has been stopped by user !");
148 fluka->SetStopEvent(kFALSE);
149 nomore = 1;
150 SOURCM.lsouit = false;
151 return;
b496f27c 152 }
153
b9d0a01d 154 //Get some info about the particle and print it
ce60a136 155 //
156 //pdg code
157 Int_t pdg = particle->GetPdgCode();
b9d0a01d 158 TVector3 polarisation;
159 particle->GetPolarisation(polarisation);
2bc4c610 160 if (debug) {
161 cout << "\t* Particle " << itrack << " retrieved..." << endl;
162 cout << "\t\t+ Name = " << particle->GetName() << endl;
163 cout << "\t\t+ PDG/Fluka code = " << pdg
4aba9d66 164 << " / " << fluka->IdFromPDG(pdg) << endl;
2bc4c610 165 cout << "\t\t+ P = ("
4aba9d66 166 << particle->Px() << " , "
167 << particle->Py() << " , "
168 << particle->Pz() << " ) --> "
169 << particle->P() << " GeV "
170 << particle->Energy() << " GeV "
171 << particle->GetMass() << " GeV " << endl;
2bc4c610 172 }
81f1d030 173 /* Npflka is the stack counter: of course any time source is called it
b9d0a01d 174 * must be =0
175 */
cfd35035 176 /* Cosines (tx,ty,tz)*/
177 Double_t cosx = particle->Px()/particle->P();
178 Double_t cosy = particle->Py()/particle->P();
160d6d40 179 Double_t cosxy = cosx * cosx + cosy * cosy;
180 Double_t cosz;
181
182 if (cosxy < 1.) {
183 cosz = TMath::Sqrt(oneone - cosxy);
184 } else {
185 cosx /= TMath::Sqrt(cosxy);
186 cosy /= TMath::Sqrt(cosxy);
187 cosz = 0.;
188 }
189
190
191
cfd35035 192 if (particle->Pz() < 0.) cosz = -cosz;
ce60a136 193
cfd35035 194 if (pdg != 50000050 && pdg != 50000051) {
4aba9d66 195 FLKSTK.npflka++;
196 Int_t ifl = fluka-> IdFromPDG(pdg);
7c33fcd7 197 Int_t ionid;
198
199 if (ifl == -2) {
200 Int_t ia = TFlukaIon::GetA(pdg);
201 Int_t iz = TFlukaIon::GetZ(pdg);
5fc9f5d5 202 Int_t is = TFlukaIon::GetIsomerNumber(pdg);
5fc9f5d5 203
204 if (is == 0) {
205 IOIOCM.iproa = ia;
206 IOIOCM.iproz = iz;
207 BEAMCM.ijhion = iz * 1000 + ia;
208 BEAMCM.ijhion = 100 * BEAMCM.ijhion + 30;
209 ionid = BEAMCM.ijhion;
210 dcdion(ionid);
211 setion(ionid);
212 FLKSTK.iloflk[FLKSTK.npflka] = ionid;
213 FLKSTK.lraddc[FLKSTK.npflka] = 0;
214 } else {
215 BEAMCM.lrdbea = 1;
216 stisbm(ia, iz, is);
217 BEAMCM.ijhion = iz * 1000 + ia;
218 BEAMCM.ijhion = 100 * BEAMCM.ijhion + 30;
219 ionid = BEAMCM.ijhion;
220 dcdion(ionid);
221 setion(ionid);
222 FLKSTK.iloflk[FLKSTK.npflka] = ionid;
223 FLKSTK.lraddc[FLKSTK.npflka] = 1;
224 }
7c33fcd7 225 } else {
226 FLKSTK.iloflk[FLKSTK.npflka] = ifl;
5fc9f5d5 227 FLKSTK.lraddc[FLKSTK.npflka] = 0;
7c33fcd7 228 ionid = ifl;
229 }
230
4aba9d66 231 /* Wtflk is the weight of the particle*/
232 FLKSTK.wtflk[FLKSTK.npflka] = oneone;
233 SUMCOU.weipri += FLKSTK.wtflk[FLKSTK.npflka];
234
235 FLKSTK.loflk[FLKSTK.npflka] = 1;
236
237 /* User dependent flag:*/
238 FLKSTK.louse[FLKSTK.npflka] = 0;
239
240 /* User dependent spare variables:*/
241 Int_t ispr = 0;
242 for (ispr = 0; ispr < mkbmx1; ispr++)
243 FLKSTK.sparek[FLKSTK.npflka][ispr] = zerzer;
244
245 /* User dependent spare flags:*/
246 for (ispr = 0; ispr < mkbmx2; ispr++)
247 FLKSTK.ispark[FLKSTK.npflka][ispr] = 0;
248
249 /* Save the track number of the stack particle:*/
250 FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1] = itrack;
251 FLKSTK.nparma++;
252 FLKSTK.numpar[FLKSTK.npflka] = FLKSTK.nparma;
253 FLKSTK.nevent[FLKSTK.npflka] = 0;
254 FLKSTK.dfnear[FLKSTK.npflka] = +zerzer;
255
256 /* Particle age (s)*/
257 FLKSTK.agestk[FLKSTK.npflka] = +zerzer;
258 FLKSTK.cmpath[FLKSTK.npflka] = +zerzer;
259 FLKSTK.aknshr[FLKSTK.npflka] = -twotwo;
260
261 /* Group number for "low" energy neutrons, set to 0 anyway*/
262 FLKSTK.igroup[FLKSTK.npflka] = 0;
263
264 /* Kinetic energy */
265 Double_t p = particle->P();
7c33fcd7 266// Double_t mass = PAPROP.am[ifl + 6];
267 Double_t mass = PAPROP.am[ionid + 6];
4aba9d66 268 FLKSTK.tkeflk[FLKSTK.npflka] = TMath::Sqrt( p * p + mass * mass) - mass;
269 /* Particle momentum*/
270 FLKSTK.pmoflk [FLKSTK.npflka] = p;
271
272 FLKSTK.txflk [FLKSTK.npflka] = cosx;
273 FLKSTK.tyflk [FLKSTK.npflka] = cosy;
274 FLKSTK.tzflk [FLKSTK.npflka] = cosz;
b9d0a01d 275
4aba9d66 276 /* Polarization cosines:*/
277 if (polarisation.Mag()) {
278 Double_t cospolx = polarisation.Px() / polarisation.Mag();
279 Double_t cospoly = polarisation.Py() / polarisation.Mag();
280 Double_t cospolz = sqrt(oneone - cospolx * cospolx - cospoly * cospoly);
281 FLKSTK.txpol [FLKSTK.npflka] = cospolx;
282 FLKSTK.typol [FLKSTK.npflka] = cospoly;
283 FLKSTK.tzpol [FLKSTK.npflka] = cospolz;
284 }
285 else {
286 FLKSTK.txpol [FLKSTK.npflka] = -twotwo;
287 FLKSTK.typol [FLKSTK.npflka] = +zerzer;
288 FLKSTK.tzpol [FLKSTK.npflka] = +zerzer;
289 }
290
291 /* Particle coordinates*/
292 // Vertext coordinates;
293 FLKSTK.xflk [FLKSTK.npflka] = particle->Vx();
294 FLKSTK.yflk [FLKSTK.npflka] = particle->Vy();
295 FLKSTK.zflk [FLKSTK.npflka] = particle->Vz();
296
297 /* Calculate the total kinetic energy of the primaries: don't change*/
298 Int_t st_ilo = FLKSTK.iloflk[FLKSTK.npflka];
299 if ( st_ilo != 0 )
300 SOURCM.tkesum +=
301 ((FLKSTK.tkeflk[FLKSTK.npflka] + PAPROP.amdisc[st_ilo+6])
302 * FLKSTK.wtflk[FLKSTK.npflka]);
303 else
304 SOURCM.tkesum += (FLKSTK.tkeflk[FLKSTK.npflka] * FLKSTK.wtflk[FLKSTK.npflka]);
305
306 /* Here we ask for the region number of the hitting point.
307 * NRGFLK (LFLKSTK) = ...
308 * The following line makes the starting region search much more
309 * robust if particles are starting very close to a boundary:
310 */
311 geocrs( FLKSTK.txflk[FLKSTK.npflka],
312 FLKSTK.tyflk[FLKSTK.npflka],
313 FLKSTK.tzflk[FLKSTK.npflka] );
b9d0a01d 314
4aba9d66 315 Int_t idisc;
316
317 georeg ( FLKSTK.xflk[FLKSTK.npflka],
318 FLKSTK.yflk[FLKSTK.npflka],
319 FLKSTK.zflk[FLKSTK.npflka],
320 FLKSTK.nrgflk[FLKSTK.npflka],
321 idisc);//<-- dummy return variable not used
322 /* Do not change these cards:*/
323 Int_t igeohsm1 = 1;
324 Int_t igeohsm2 = -11;
325 geohsm ( FLKSTK.nhspnt[FLKSTK.npflka], igeohsm1, igeohsm2, LTCLCM.mlattc );
326 FLKSTK.nlattc[FLKSTK.npflka] = LTCLCM.mlattc;
327 soevsv();
cfd35035 328 } else {
b496f27c 329 //
4aba9d66 330 // Next particle is optical photon
331 //
332 OPPHST.lstopp++;
333 OPPHST.donear [OPPHST.lstopp - 1] = 0.;
334
335 OPPHST.xoptph [OPPHST.lstopp - 1] = particle->Vx();
336 OPPHST.yoptph [OPPHST.lstopp - 1] = particle->Vy();
337 OPPHST.zoptph [OPPHST.lstopp - 1] = particle->Vz();
338
339 OPPHST.txopph [OPPHST.lstopp - 1] = cosx;
340 OPPHST.tyopph [OPPHST.lstopp - 1] = cosy;
341 OPPHST.tzopph [OPPHST.lstopp - 1] = cosz;
342
343
344 if (polarisation.Mag()) {
345 Double_t cospolx = polarisation.Px() / polarisation.Mag();
346 Double_t cospoly = polarisation.Py() / polarisation.Mag();
347 Double_t cospolz = sqrt(oneone - cospolx * cospolx - cospoly * cospoly);
348 OPPHST.txpopp [OPPHST.lstopp - 1] = cospolx;
349 OPPHST.typopp [OPPHST.lstopp - 1] = cospoly;
350 OPPHST.tzpopp [OPPHST.lstopp - 1] = cospolz;
351 }
352 else {
353 OPPHST.txpopp [OPPHST.lstopp - 1] = -twotwo;
354 OPPHST.typopp [OPPHST.lstopp - 1] = +zerzer;
355 OPPHST.tzpopp [OPPHST.lstopp - 1] = +zerzer;
356 }
357
358 geocrs( OPPHST.txopph[OPPHST.lstopp - 1],
359 OPPHST.tyopph[OPPHST.lstopp - 1],
360 OPPHST.tzopph[OPPHST.lstopp - 1] );
361
362 Int_t idisc;
363
364 georeg ( OPPHST.xoptph[OPPHST.lstopp - 1],
365 OPPHST.yoptph[OPPHST.lstopp - 1],
366 OPPHST.zoptph[OPPHST.lstopp - 1],
367 OPPHST.nregop[OPPHST.lstopp - 1],
368 idisc);//<-- dummy return variable not used
369
370 OPPHST.wtopph [OPPHST.lstopp - 1] = particle->GetWeight();
371 OPPHST.poptph [OPPHST.lstopp - 1] = particle->P();
372 OPPHST.agopph [OPPHST.lstopp - 1] = particle->T();
373 OPPHST.cmpopp [OPPHST.lstopp - 1] = +zerzer;
374 OPPHST.loopph [OPPHST.lstopp - 1] = 0;
375 OPPHST.louopp [OPPHST.lstopp - 1] = itrack;
376 OPPHST.nlatop [OPPHST.lstopp - 1] = LTCLCM.mlattc;
377 }
378
ce60a136 379//
380// Pre-track actions at for primary tracks
381//
382 if (particleIsPrimary) {
5cc3d37e 383 fluka->SetCaller(kSODRAW);
4aba9d66 384 TVirtualMCApplication::Instance()->BeginPrimary();
385 TVirtualMCApplication::Instance()->PreTrack();
ce60a136 386 }
ce60a136 387//
2bc4c610 388 if (debug) cout << "<== source(" << nomore << ")" << endl;
b9d0a01d 389 }
390}