User stepping methods added (E. Futo)
[u/mrichter/AliRoot.git] / TFluka / TFluka.cxx
CommitLineData
03ca248b 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
17$Log$
fa3d1cc7 18Revision 1.6 2002/11/21 18:40:06 iglez2
19Media handling added
20
27b2f7fe 21Revision 1.5 2002/11/07 17:59:10 iglez2
22Included the geometry through geant4_vmc/FLUGG
23
bf3aa28e 24Revision 1.4 2002/11/04 16:00:46 iglez2
25The conversion between ID and PDG now uses Fluka routines and arrays which is more consistent.
26
f9cb2fec 27Revision 1.3 2002/10/22 15:12:14 alibrary
28Introducing Riostream.h
29
eae0fe66 30Revision 1.2 2002/10/14 14:57:40 hristov
31Merging the VirtualMC branch to the main development branch (HEAD)
32
b9d0a01d 33Revision 1.1.2.8 2002/10/08 16:33:17 iglez2
34LSOUIT is set to true before the second call to flukam.
35
36Revision 1.1.2.7 2002/10/08 09:30:37 iglez2
37Solved stupid missing ;
38
39Revision 1.1.2.6 2002/10/07 13:40:22 iglez2
40First implementations of the PDG <--> Fluka Id conversion routines
41
42Revision 1.1.2.5 2002/09/26 16:26:03 iglez2
43Added verbosity
44Call to gAlice->Generator()->Generate()
45
46Revision 1.1.2.4 2002/09/26 13:22:23 iglez2
47Naive implementation of ProcessRun and ProcessEvent
48Opening/Closing of input file (fInputFileName) with FORTRAN unit 5 before/after the first call to flukam inside Init()
49
50Revision 1.1.2.3 2002/09/20 15:35:51 iglez2
51Modification of LFDRTR. Value is passed to FLUKA !!!
52
53Revision 1.1.2.2 2002/09/18 14:34:44 iglez2
54Revised version with all pure virtual methods implemented
55
56Revision 1.1.2.1 2002/07/24 08:49:41 alibrary
57Adding TFluka to VirtualMC
58
59Revision 1.1 2002/07/05 13:10:07 morsch
60First commit of Fluka interface.
61
03ca248b 62*/
63
eae0fe66 64#include <Riostream.h>
b9d0a01d 65
03ca248b 66#include "TFluka.h"
b9d0a01d 67#include "TCallf77.h" //For the fortran calls
68#include "Fdblprc.h" //(DBLPRC) fluka common
b9d0a01d 69#include "Fepisor.h" //(EPISOR) fluka common
fa3d1cc7 70#include "Ffinuc.h" //(FINUC) fluka common
71#include "Fiounit.h" //(IOUNIT) fluka common
72#include "Fpaprop.h" //(PAPROP) fluka common
f9cb2fec 73#include "Fpart.h" //(PART) fluka common
fa3d1cc7 74#include "Ftrackr.h" //(TRACKR) fluka common
75#include "Ffheavy.h" //(FHEAVY) fluka common
b9d0a01d 76
fa3d1cc7 77#include "TVirtualMC.h"
bf3aa28e 78#include "TG4GeometryManager.h" //For the geometry management
79#include "TG4DetConstruction.h" //For the detector construction
80
81#include "FGeometryInit.hh"
fa3d1cc7 82#include "TLorentzVector.h"
bf3aa28e 83
b9d0a01d 84// Fluka methods that may be needed.
85#ifndef WIN32
86# define flukam flukam_
87# define fluka_openinp fluka_openinp_
88# define fluka_closeinp fluka_closeinp_
f9cb2fec 89# define mcihad mcihad_
90# define mpdgha mpdgha_
b9d0a01d 91#else
92# define flukam FLUKAM
93# define fluka_openinp FLUKA_OPENINP
94# define fluka_closeinp FLUKA_CLOSEINP
f9cb2fec 95# define mcihad MCIHAD
96# define mpdgha MPDGHA
b9d0a01d 97#endif
98
99extern "C"
100{
101 //
102 // Prototypes for FLUKA functions
103 //
104 void type_of_call flukam(const int&);
105 void type_of_call fluka_openinp(const int&, DEFCHARA);
106 void type_of_call fluka_closeinp(const int&);
f9cb2fec 107 int type_of_call mcihad(const int&);
108 int type_of_call mpdgha(const int&);
b9d0a01d 109}
110
111//
112// Class implementation for ROOT
113//
03ca248b 114ClassImp(TFluka)
b9d0a01d 115
116//
bf3aa28e 117//----------------------------------------------------------------------------
118// TFluka constructors and destructors.
b9d0a01d 119//____________________________________________________________________________
120TFluka::TFluka()
121 :TVirtualMC(),
122 fVerbosityLevel(0),
bf3aa28e 123 fInputFileName(""),
27b2f7fe 124 fDetector(0),
125 fCurrentFlukaRegion(-1)
b9d0a01d 126{
127 //
128 // Default constructor
129 //
130}
131
b9d0a01d 132TFluka::TFluka(const char *title, Int_t verbosity)
133 :TVirtualMC("TFluka",title),
134 fVerbosityLevel(verbosity),
bf3aa28e 135 fInputFileName(""),
27b2f7fe 136 fDetector(0),
137 fCurrentFlukaRegion(-1)
b9d0a01d 138{
139 if (fVerbosityLevel >=3)
140 cout << "==> TFluka::TFluka(" << title << ") constructor called." << endl;
141
bf3aa28e 142
143 // create geometry manager
144 if (fVerbosityLevel >=2)
145 cout << "\t* Creating G4 Geometry manager..." << endl;
146 fGeometryManager = new TG4GeometryManager();
147 if (fVerbosityLevel >=2)
148 cout << "\t* Creating G4 Detector..." << endl;
149 fDetector = new TG4DetConstruction();
150 FGeometryInit* geominit = FGeometryInit::GetInstance();
151 if (geominit)
152 geominit->setDetConstruction(fDetector);
153 else {
154 cerr << "ERROR: Could not create FGeometryInit!" << endl;
155 cerr << " Exiting!!!" << endl;
156 abort();
157 }
158
b9d0a01d 159 if (fVerbosityLevel >=3)
160 cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
161}
162
bf3aa28e 163TFluka::~TFluka() {
164 if (fVerbosityLevel >=3)
165 cout << "==> TFluka::~TFluka() destructor called." << endl;
166
167 delete fGeometryManager;
168
169 if (fVerbosityLevel >=3)
170 cout << "<== TFluka::~TFluka() destructor called." << endl;
171}
172
173//
174//_____________________________________________________________________________
175// TFluka control methods
b9d0a01d 176//____________________________________________________________________________
177void TFluka::Init() {
178 if (fVerbosityLevel >=3)
179 cout << "==> TFluka::Init() called." << endl;
180
181 if (fVerbosityLevel >=2)
182 cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
183 << ") in fluka..." << endl;
184 GLOBAL.lfdrtr = true;
185
186 if (fVerbosityLevel >=2)
187 cout << "\t* Opening file " << fInputFileName << endl;
188 const char* fname = fInputFileName;
189 fluka_openinp(lunin, PASSCHARA(fname));
190
191 if (fVerbosityLevel >=2)
192 cout << "\t* Calling flukam..." << endl;
bf3aa28e 193 flukam(1);
b9d0a01d 194
195 if (fVerbosityLevel >=2)
196 cout << "\t* Closing file " << fInputFileName << endl;
197 fluka_closeinp(lunin);
198
199 if (fVerbosityLevel >=3)
200 cout << "<== TFluka::Init() called." << endl;
fa3d1cc7 201
b9d0a01d 202}
203
bf3aa28e 204void TFluka::FinishGeometry() {
205 if (fVerbosityLevel >=3)
206 cout << "==> TFluka::FinishGeometry() called." << endl;
207
208 fGeometryManager->Ggclos();
209
27b2f7fe 210 FGeometryInit* flugg = FGeometryInit::GetInstance();
211 map<TString, Int_t, less<TString> >::iterator i;
212 for (fVolumeMediaMap.begin(); i != fVolumeMediaMap.end(); i++) {
213 TString volName = (*i).first;
214 Int_t media = (*i).second;
215 Int_t region = flugg->GetRegionFromName(volName);
216 fMediaByRegion[region] = media;
217 }
218
bf3aa28e 219 if (fVerbosityLevel >=3)
220 cout << "<== TFluka::FinishGeometry() called." << endl;
221}
222
223void TFluka::BuildPhysics() {
224 if (fVerbosityLevel >=3)
225 cout << "==> TFluka::BuildPhysics() called." << endl;
226
227
228 if (fVerbosityLevel >=3)
229 cout << "<== TFluka::BuildPhysics() called." << endl;
230}
231
b9d0a01d 232void TFluka::ProcessEvent() {
233 if (fVerbosityLevel >=3)
234 cout << "==> TFluka::ProcessEvent() called." << endl;
235
236 if (fVerbosityLevel >=3)
237 cout << "<== TFluka::ProcessEvent() called." << endl;
238}
239
bf3aa28e 240
b9d0a01d 241void TFluka::ProcessRun(Int_t nevent) {
242 if (fVerbosityLevel >=3)
243 cout << "==> TFluka::ProcessRun(" << nevent << ") called."
244 << endl;
245
246 if (fVerbosityLevel >=2) {
247 cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
248 cout << "\t* Calling flukam again..." << endl;
249 }
250 fApplication->GeneratePrimaries();
251 EPISOR.lsouit = true;
bf3aa28e 252 flukam(1);
b9d0a01d 253
254 if (fVerbosityLevel >=3)
255 cout << "<== TFluka::ProcessRun(" << nevent << ") called."
256 << endl;
257}
258
bf3aa28e 259//_____________________________________________________________________________
260// methods for building/management of geometry
261//____________________________________________________________________________
262// functions from GCONS
263void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
264 Float_t &dens, Float_t &radl, Float_t &absl,
265 Float_t* ubuf, Int_t& nbuf) {
266//
267 fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
268}
269
270void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
271 Double_t &dens, Double_t &radl, Double_t &absl,
272 Double_t* ubuf, Int_t& nbuf) {
273//
274 fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
275}
276
277// detector composition
278void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
279 Double_t z, Double_t dens, Double_t radl, Double_t absl,
280 Float_t* buf, Int_t nwbuf) {
281//
282 fGeometryManager
283 ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
284}
285void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
286 Double_t z, Double_t dens, Double_t radl, Double_t absl,
287 Double_t* buf, Int_t nwbuf) {
288//
289 fGeometryManager
290 ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
291}
292
293void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
294 Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
295//
296 fGeometryManager
297 ->Mixture(kmat, name, a, z, dens, nlmat, wmat);
298}
299void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
300 Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
301//
302 fGeometryManager
303 ->Mixture(kmat, name, a, z, dens, nlmat, wmat);
304}
305
306void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
307 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
308 Double_t stemax, Double_t deemax, Double_t epsil,
309 Double_t stmin, Float_t* ubuf, Int_t nbuf) {
310 //
311 fGeometryManager
312 ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
313 epsil, stmin, ubuf, nbuf);
314}
315void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
316 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
317 Double_t stemax, Double_t deemax, Double_t epsil,
318 Double_t stmin, Double_t* ubuf, Int_t nbuf) {
319 //
320 fGeometryManager
321 ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
322 epsil, stmin, ubuf, nbuf);
323}
324
325void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
326 Double_t thetaY, Double_t phiY, Double_t thetaZ,
327 Double_t phiZ) {
328//
329 fGeometryManager
330 ->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
331}
332
333void TFluka::Gstpar(Int_t itmed, const char *param, Double_t parval) {
334//
335 fGeometryManager->Gstpar(itmed, param, parval);
336}
337
338// functions from GGEOM
339Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
340 Float_t *upar, Int_t np) {
341//
27b2f7fe 342 fVolumeMediaMap[TString(name)] = nmed;
bf3aa28e 343 return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
344}
345Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
346 Double_t *upar, Int_t np) {
347//
348 return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
349}
350
351void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
352 Int_t iaxis) {
353//
354 fGeometryManager->Gsdvn(name, mother, ndiv, iaxis);
355}
356
357void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
358 Int_t iaxis, Double_t c0i, Int_t numed) {
359//
360 fGeometryManager->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
361}
362
363void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
364 Int_t iaxis, Int_t numed, Int_t ndvmx) {
365//
366 fGeometryManager->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
367}
368
369void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
370 Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
371//
372 fGeometryManager->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
373}
374
375void TFluka::Gsord(const char *name, Int_t iax) {
376//
377 fGeometryManager->Gsord(name, iax);
378}
379
380void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
381 Double_t x, Double_t y, Double_t z, Int_t irot,
382 const char *konly) {
383//
384 fGeometryManager->Gspos(name, nr, mother, x, y, z, irot, konly);
385}
386
387void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
388 Double_t x, Double_t y, Double_t z, Int_t irot,
389 const char *konly, Float_t *upar, Int_t np) {
390 //
391 fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
392}
393void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
394 Double_t x, Double_t y, Double_t z, Int_t irot,
395 const char *konly, Double_t *upar, Int_t np) {
396 //
397 fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
398}
399
400void TFluka::Gsbool(const char* onlyVolName, const char* manyVolName) {
401//
402 fGeometryManager->Gsbool(onlyVolName, manyVolName);
403}
404
405void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t *ppckov,
406 Float_t *absco, Float_t *effic, Float_t *rindex) {
407//
408 fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
409}
410void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
411 Double_t *absco, Double_t *effic, Double_t *rindex) {
412//
413 fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
414}
415
416// Euclid
417void TFluka::WriteEuclid(const char* fileName, const char* topVol,
418 Int_t number, Int_t nlevel) {
419//
420 fGeometryManager->WriteEuclid(fileName, topVol, number, nlevel);
421}
422
423
424
27b2f7fe 425//_____________________________________________________________________________
426// methods needed by the stepping
427//____________________________________________________________________________
428Int_t TFluka::GetMedium() const {
429 return fMediaByRegion[fCurrentFlukaRegion];
430}
bf3aa28e 431
432
433
434//____________________________________________________________________________
435// ID <--> PDG transformations
b9d0a01d 436//_____________________________________________________________________________
437Int_t TFluka::IdFromPDG(Int_t pdg) const
438{
439 //
f9cb2fec 440 // Return Fluka code from PDG and pseudo ENDF code
441
442 // MCIHAD() goes from pdg to fluka internal.
443 Int_t intfluka = mcihad(pdg);
444 // KPTOIP array goes from internal to official
445 return GetFlukaKPTOIP(intfluka);
b9d0a01d 446}
447
b9d0a01d 448Int_t TFluka::PDGFromId(Int_t id) const
449{
450 //
f9cb2fec 451 // Return PDG code and pseudo ENDF code from Fluka code
b9d0a01d 452
f9cb2fec 453 //IPTOKP array goes from official to internal
454 Int_t intfluka = GetFlukaIPTOKP(id);
455 //MPKDHA() goes from internal to PDG
456 return mpdgha(intfluka);
457
b9d0a01d 458}
fa3d1cc7 459
460
461
462//_____________________________________________________________________________
463// methods for step management
464//____________________________________________________________________________
465//
466// dynamic properties
467//
468void TFluka::TrackPosition(TLorentzVector& position) const
469{
470// Return the current position in the master reference frame of the
471// track being transported
472// TRACKR.atrack = age of the particle
473// TRACKR.xtrack = x-position of the last point
474// TRACKR.ytrack = y-position of the last point
475// TRACKR.ztrack = z-position of the last point
476 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
477 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
478 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
479 position.SetT(TRACKR.atrack);
480}
481
482void TFluka::TrackMomentum(TLorentzVector& momentum) const
483{
484// Return the direction and the momentum (GeV/c) of the track
485// currently being transported
486// TRACKR.ptrack = momentum of the particle (not always defined, if
487// < 0 must be obtained from etrack)
488// TRACKR.cx,y,ztrck = direction cosines of the current particle
489// TRACKR.etrack = total energy of the particle
490// TRACKR.jtrack = identity number of the particle
491// PAPROP.am[TRACKR.jtrack] = particle mass in gev
492 if (TRACKR.ptrack >= 0) {
493 momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
494 momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
495 momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
496 momentum.SetE(TRACKR.etrack);
497 return;
498 }
499 else {
500 Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack]*PAPROP.am[TRACKR.jtrack]);
501 momentum.SetPx(p*TRACKR.cxtrck);
502 momentum.SetPy(p*TRACKR.cytrck);
503 momentum.SetPz(p*TRACKR.cztrck);
504 momentum.SetE(TRACKR.etrack);
505 return;
506 }
507}
508
509Double_t TFluka::TrackStep() const
510{
511// Return the length in centimeters of the current step
512// TRACKR.ctrack = total curved path
513 return TRACKR.ctrack;
514}
515
516Double_t TFluka::TrackLength() const
517{
518// It is wrong
519// should be the sum of all steps starting from the beginning of the track
520// for the time being returns only the length in centimeters of the current step
521 return TRACKR.ctrack;
522}
523
524Double_t TFluka::TrackTime() const
525{
526// Return the current time of flight of the track being transported
527// TRACKR.atrack = age of the particle
528 return TRACKR.atrack;
529}
530
531Double_t TFluka::Edep() const
532{
533// Energy deposition
534// if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
535// -->local energy deposition (the value and the point are not recorded in TRACKR)
536// but in the variable "rull" of the procedure "endraw.cxx"
537// if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
538// -->no energy loss along the track
539// if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
540// -->energy loss distributed along the track
541// TRACKR.dtrack = energy deposition of the jth deposition even
542 if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
543 return fRull;
544 else {
545 Double_t sum = 0;
546 for ( Int_t j=0;j<TRACKR.mtrack;j++) {
547 sum +=TRACKR.dtrack[j];
548 }
549 return sum;
550 }
551}
552
553Int_t TFluka::TrackPid() const
554{
555// Return the id of the particle transported
556// TRACKR.jtrack = identity number of the particle
557 return PDGFromId(TRACKR.jtrack);
558}
559
560Double_t TFluka::TrackCharge() const
561{
562// Return charge of the track currently transported
563// PAPROP.ichrge = electric charge of the particle
564 return PAPROP.ichrge[TRACKR.jtrack+6];
565}
566
567Double_t TFluka::TrackMass() const
568{
569// PAPROP.am = particle mass in GeV
570 return PAPROP.am[TRACKR.jtrack+6];
571}
572
573Double_t TFluka::Etot() const
574{
575// TRACKR.etrack = total energy of the particle
576 return TRACKR.etrack;
577}
578
579//
580// track status
581//
582Bool_t TFluka::IsNewTrack() const
583{
584// ???????????????,
585// True if the track is not at the boundary of the current volume
586// Not true in some cases in bxdraw - to be solved
587 return 1;
588}
589
590Bool_t TFluka::IsTrackInside() const
591{
592// True if the track is not at the boundary of the current volume
593// In Fluka a step is always inside one kind of material
594// If the step would go behind the region of one material,
595// it will be shortened to reach only the boundary.
596// Therefore IsTrackInside() is always true.
597// Not true in some cases in bxdraw - to be solved
598 return 1;
599}
600
601Bool_t TFluka::IsTrackEntering() const
602{
603// True if this is the first step of the track in the current volume
604// Boundary- (X) crossing
605// Icode = 19: boundary crossing - call from Kaskad
606// Icode = 29: boundary crossing - call from Emfsco
607// Icode = 39: boundary crossing - call from Kasneu
608// Icode = 49: boundary crossing - call from Kashea
609// Icode = 59: boundary crossing - call from Kasoph
610 if (fIcode == 19 ||
611 fIcode == 29 ||
612 fIcode == 39 ||
613 fIcode == 49 ||
614 fIcode == 59) return 1;
615 else return 0;
616}
617
618Bool_t TFluka::IsTrackExiting() const
619{
620// True if this is the last step of the track in the current volume
621// Boundary- (X) crossing
622// Icode = 19: boundary crossing - call from Kaskad
623// Icode = 29: boundary crossing - call from Emfsco
624// Icode = 39: boundary crossing - call from Kasneu
625// Icode = 49: boundary crossing - call from Kashea
626// Icode = 59: boundary crossing - call from Kasoph
627 if (fIcode == 19 ||
628 fIcode == 29 ||
629 fIcode == 39 ||
630 fIcode == 49 ||
631 fIcode == 59) return 1;
632 else return 0;
633}
634
635Bool_t TFluka::IsTrackOut() const
636{
637// True if the track is out of the setup
638// means escape
639// Icode = 14: escape - call from Kaskad
640// Icode = 23: escape - call from Emfsco
641// Icode = 32: escape - call from Kasneu
642// Icode = 40: escape - call from Kashea
643// Icode = 51: escape - call from Kasoph
644 if (fIcode == 14 ||
645 fIcode == 23 ||
646 fIcode == 32 ||
647 fIcode == 40 ||
648 fIcode == 51) return 1;
649 else return 0;
650}
651
652Bool_t TFluka::IsTrackDisappeared() const
653{
654// means all inelastic interactions and decays
655// fIcode from usdraw
656 if (fIcode == 101 || // inelastic interaction
657 fIcode == 102 || // particle decay
658 fIcode == 214 || // in-flight annihilation
659 fIcode == 215 || // annihilation at rest
660 fIcode == 217 || // pair production
661 fIcode == 221) return 1;
662 else return 0;
663}
664
665Bool_t TFluka::IsTrackStop() const
666{
667// True if the track energy has fallen below the threshold
668// means stopped by signal or below energy threshold
669// Icode = 12: stopping particle - call from Kaskad
670// Icode = 15: time kill - call from Kaskad
671// Icode = 21: below threshold, iarg=1 - call from Emfsco
672// Icode = 22: below threshold, iarg=2 - call from Emfsco
673// Icode = 24: time kill - call from Emfsco
674// Icode = 31: below threshold - call from Kasneu
675// Icode = 33: time kill - call from Kasneu
676// Icode = 41: time kill - call from Kashea
677// Icode = 52: time kill - call from Kasoph
678 if (fIcode == 12 ||
679 fIcode == 15 ||
680 fIcode == 21 ||
681 fIcode == 22 ||
682 fIcode == 24 ||
683 fIcode == 31 ||
684 fIcode == 33 ||
685 fIcode == 41 ||
686 fIcode == 52) return 1;
687 else return 0;
688}
689
690Bool_t TFluka::IsTrackAlive() const
691{
692// means not disappeared or not out
693 if (IsTrackDisappeared() || IsTrackOut() ) return 0;
694 else return 1;
695}
696
697//
698// secondaries
699//
700
701Int_t TFluka::NSecondaries() const
702// Number of secondary particles generated in the current step
703// fIcode = 100 = elastic interaction
704// fIcode = 101 = inelastic interaction
705// fIcode = 102 = particle decay
706// fIcode = 103 = delta ray
707// fIcode = 104 = pair production
708// fIcode = 105 = bremsstrahlung
709{
710 if (fIcode >= 100 && fIcode <= 105)
711 return FINUC.np + FHEAVY.npheav;
712 else
713 return -1;
714}
715
716void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
717 TLorentzVector& position, TLorentzVector& momentum)
718//
719// fIcode = 100 = elastic interaction
720// fIcode = 101 = inelastic interaction
721// fIcode = 102 = particle decay
722// fIcode = 103 = delta ray
723// fIcode = 104 = pair production
724// fIcode = 105 = bremsstrahlung
725{
726 if (fIcode >= 100 && fIcode <= 105) {
727 if (isec >= 0 && isec < FINUC.np) {
728 particleId = PDGFromId(FINUC.kpart[isec]);
729 position.SetX(fXsco);
730 position.SetY(fYsco);
731 position.SetZ(fZsco);
732 position.SetT(FINUC.agesec[isec]);
733 momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
734 momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
735 momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
736 momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
737 }
738 if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
739 Int_t jsec = isec - FINUC.np;
740 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
741 position.SetX(fXsco);
742 position.SetY(fYsco);
743 position.SetZ(fZsco);
744 position.SetT(FHEAVY.agheav[jsec]);
745 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
746 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
747 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
748 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
749 momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
750 else if (FHEAVY.tkheav[jsec] > 6)
751 momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
752 }
753 }
754}
755
756//TMCProcess ProdProcess(Int_t isec) const
757// Name of the process that has produced the secondary particles
758// in the current step
759//{
760// will come from FINUC when called from USDRAW
761//}
762
763//Int_t StepProcesses(TArrayI &proc) const
764// Return processes active in the current step
765//{
766//ck = total energy of the particl ????????????????
767//}
768
769
770// ===============================================================
771void TFluka::FutoTest()
772{
773 Int_t icode, mreg, newreg, particleId;
774// Int_t medium;
775 Double_t rull, xsco, ysco, zsco;
776 TLorentzVector position, momentum;
777 icode = GetIcode();
778 if (icode == 0) {
779 cout << " icode=" << icode << endl;
780 /*
781 cout << "TLorentzVector positionX=" << position.X()
782 << "positionY=" << position.Y()
783 << "positionZ=" << position.Z()
784 << "timeT=" << position.T() << endl;
785 cout << "TLorentzVector momentumX=" << momentum.X()
786 << "momentumY=" << momentum.Y()
787 << "momentumZ=" << momentum.Z()
788 << "energyE=" << momentum.E() << endl;
789 cout << "TrackPid=" << TrackPid() << endl;
790 */
791 }
792
793 else if (icode > 0 && icode <= 5) {
794// mgdraw
795 mreg = GetMreg();
796// medium = GetMedium();
797 cout << " icode=" << icode
798 << " mreg=" << mreg
799// << " medium=" << medium
800 << endl;
801 TrackPosition(position);
802 TrackMomentum(momentum);
803 cout << "TLorentzVector positionX=" << position.X()
804 << "positionY=" << position.Y()
805 << "positionZ=" << position.Z()
806 << "timeT=" << position.T() << endl;
807 cout << "TLorentzVector momentumX=" << momentum.X()
808 << "momentumY=" << momentum.Y()
809 << "momentumZ=" << momentum.Z()
810 << "energyE=" << momentum.E() << endl;
811 cout << "TrackStep=" << TrackStep() << endl;
812 cout << "TrackLength=" << TrackLength() << endl;
813 cout << "TrackTime=" << TrackTime() << endl;
814 cout << "Edep=" << Edep() << endl;
815 cout << "TrackPid=" << TrackPid() << endl;
816 cout << "TrackCharge=" << TrackCharge() << endl;
817 cout << "TrackMass=" << TrackMass() << endl;
818 cout << "Etot=" << Etot() << endl;
819 cout << "IsNewTrack=" << IsNewTrack() << endl;
820 cout << "IsTrackInside=" << IsTrackInside() << endl;
821 cout << "IsTrackEntering=" << IsTrackEntering() << endl;
822 cout << "IsTrackExiting=" << IsTrackExiting() << endl;
823 cout << "IsTrackOut=" << IsTrackOut() << endl;
824 cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
825 cout << "IsTrackAlive=" << IsTrackAlive() << endl;
826 }
827
828 else if((icode >= 10 && icode <= 15) ||
829 (icode >= 20 && icode <= 24) ||
830 (icode >= 30 && icode <= 33) ||
831 (icode >= 40 && icode <= 41) ||
832 (icode >= 50 && icode <= 52)) {
833// endraw
834 mreg = GetMreg();
835// medium = GetMedium();
836 rull = GetRull();
837 xsco = GetXsco();
838 ysco = GetYsco();
839 zsco = GetZsco();
840 cout << " icode=" << icode
841 << " mreg=" << mreg
842// << " medium=" << medium
843 << " rull=" << rull
844 << " xsco=" << xsco
845 << " ysco=" << ysco
846 << " zsco=" << zsco << endl;
847 TrackPosition(position);
848 TrackMomentum(momentum);
849 cout << "Edep=" << Edep() << endl;
850 cout << "Etot=" << Etot() << endl;
851 cout << "TrackPid=" << TrackPid() << endl;
852 cout << "TrackCharge=" << TrackCharge() << endl;
853 cout << "TrackMass=" << TrackMass() << endl;
854 cout << "IsTrackOut=" << IsTrackOut() << endl;
855 cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
856 cout << "IsTrackStop=" << IsTrackStop() << endl;
857 cout << "IsTrackAlive=" << IsTrackAlive() << endl;
858 }
859
860 else if((icode >= 100 && icode <= 105) ||
861 (icode == 208) ||
862 (icode == 210) ||
863 (icode == 212) ||
864 (icode >= 214 && icode <= 215) ||
865 (icode == 217) ||
866 (icode == 219) ||
867 (icode == 221) ||
868 (icode == 225) ||
869 (icode == 300) ||
870 (icode == 400)) {
871// usdraw
872 mreg = GetMreg();
873// medium = GetMedium();
874 xsco = GetXsco();
875 ysco = GetYsco();
876 zsco = GetZsco();
877 cout << " icode=" << icode
878 << " mreg=" << mreg
879// << " medium=" << medium
880 << " xsco=" << xsco
881 << " ysco=" << ysco
882 << " zsco=" << zsco << endl;
883 cout << "TrackPid=" << TrackPid() << endl;
884 cout << "NSecondaries=" << NSecondaries() << endl;
885 for (Int_t isec=0; isec< NSecondaries(); isec++) {
886//void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
887// TLorentzVector& position, TLorentzVector& momentum)
888 TFluka::GetSecondary(isec, particleId, position, momentum);
889 cout << "TLorentzVector positionX=" << position.X()
890 << "positionY=" << position.Y()
891 << "positionZ=" << position.Z()
892 << "timeT=" << position.T() << endl;
893 cout << "TLorentzVector momentumX=" << momentum.X()
894 << "momentumY=" << momentum.Y()
895 << "momentumZ=" << momentum.Z()
896 << "energyE=" << momentum.E() << endl;
897 cout << "TrackPid=" << particleId << endl;
898
899 }
900 }
901
902 else if((icode == 19) ||
903 (icode == 29) ||
904 (icode == 39) ||
905 (icode == 49) ||
906 (icode == 59)) {
907 mreg = GetMreg();
908// medium = GetMedium();
909 newreg = GetNewreg();
910 xsco = GetXsco();
911 ysco = GetYsco();
912 zsco = GetZsco();
913 cout << " icode=" << icode
914 << " mreg=" << mreg
915// << " medium=" << medium
916 << " newreg=" << newreg
917 << " xsco=" << xsco
918 << " ysco=" << ysco
919 << " zsco=" << zsco << endl;
920 }
921//
922// ====================================================================
923//
924
925
926
927} // end of FutoTest
928