1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.7 2002/12/06 12:21:32 morsch
19 User stepping methods added (E. Futo)
21 Revision 1.6 2002/11/21 18:40:06 iglez2
24 Revision 1.5 2002/11/07 17:59:10 iglez2
25 Included the geometry through geant4_vmc/FLUGG
27 Revision 1.4 2002/11/04 16:00:46 iglez2
28 The conversion between ID and PDG now uses Fluka routines and arrays which is more consistent.
30 Revision 1.3 2002/10/22 15:12:14 alibrary
31 Introducing Riostream.h
33 Revision 1.2 2002/10/14 14:57:40 hristov
34 Merging the VirtualMC branch to the main development branch (HEAD)
36 Revision 1.1.2.8 2002/10/08 16:33:17 iglez2
37 LSOUIT is set to true before the second call to flukam.
39 Revision 1.1.2.7 2002/10/08 09:30:37 iglez2
40 Solved stupid missing ;
42 Revision 1.1.2.6 2002/10/07 13:40:22 iglez2
43 First implementations of the PDG <--> Fluka Id conversion routines
45 Revision 1.1.2.5 2002/09/26 16:26:03 iglez2
47 Call to gAlice->Generator()->Generate()
49 Revision 1.1.2.4 2002/09/26 13:22:23 iglez2
50 Naive implementation of ProcessRun and ProcessEvent
51 Opening/Closing of input file (fInputFileName) with FORTRAN unit 5 before/after the first call to flukam inside Init()
53 Revision 1.1.2.3 2002/09/20 15:35:51 iglez2
54 Modification of LFDRTR. Value is passed to FLUKA !!!
56 Revision 1.1.2.2 2002/09/18 14:34:44 iglez2
57 Revised version with all pure virtual methods implemented
59 Revision 1.1.2.1 2002/07/24 08:49:41 alibrary
60 Adding TFluka to VirtualMC
62 Revision 1.1 2002/07/05 13:10:07 morsch
63 First commit of Fluka interface.
67 #include <Riostream.h>
69 #include "TClonesArray.h"
71 #include "TCallf77.h" //For the fortran calls
72 #include "Fdblprc.h" //(DBLPRC) fluka common
73 #include "Fepisor.h" //(EPISOR) fluka common
74 #include "Ffinuc.h" //(FINUC) fluka common
75 #include "Fiounit.h" //(IOUNIT) fluka common
76 #include "Fpaprop.h" //(PAPROP) fluka common
77 #include "Fpart.h" //(PART) fluka common
78 #include "Ftrackr.h" //(TRACKR) fluka common
79 #include "Fpaprop.h" //(PAPROP) fluka common
80 #include "Ffheavy.h" //(FHEAVY) fluka common
82 #include "TVirtualMC.h"
83 #include "TG4GeometryManager.h" //For the geometry management
84 #include "TG4DetConstruction.h" //For the detector construction
86 #include "FGeometryInit.hh"
87 #include "TLorentzVector.h"
88 #include "FlukaVolume.h"
90 // Fluka methods that may be needed.
92 # define flukam flukam_
93 # define fluka_openinp fluka_openinp_
94 # define fluka_closeinp fluka_closeinp_
95 # define mcihad mcihad_
96 # define mpdgha mpdgha_
98 # define flukam FLUKAM
99 # define fluka_openinp FLUKA_OPENINP
100 # define fluka_closeinp FLUKA_CLOSEINP
101 # define mcihad MCIHAD
102 # define mpdgha MPDGHA
108 // Prototypes for FLUKA functions
110 void type_of_call flukam(const int&);
111 void type_of_call fluka_openinp(const int&, DEFCHARA);
112 void type_of_call fluka_closeinp(const int&);
113 int type_of_call mcihad(const int&);
114 int type_of_call mpdgha(const int&);
118 // Class implementation for ROOT
123 //----------------------------------------------------------------------------
124 // TFluka constructors and destructors.
125 //____________________________________________________________________________
131 fCurrentFlukaRegion(-1)
134 // Default constructor
138 TFluka::TFluka(const char *title, Int_t verbosity)
139 :TVirtualMC("TFluka",title),
140 fVerbosityLevel(verbosity),
143 fCurrentFlukaRegion(-1)
145 if (fVerbosityLevel >=3)
146 cout << "==> TFluka::TFluka(" << title << ") constructor called." << endl;
149 // create geometry manager
150 if (fVerbosityLevel >=2)
151 cout << "\t* Creating G4 Geometry manager..." << endl;
152 fGeometryManager = new TG4GeometryManager();
153 if (fVerbosityLevel >=2)
154 cout << "\t* Creating G4 Detector..." << endl;
155 fDetector = new TG4DetConstruction();
156 FGeometryInit* geominit = FGeometryInit::GetInstance();
158 geominit->setDetConstruction(fDetector);
160 cerr << "ERROR: Could not create FGeometryInit!" << endl;
161 cerr << " Exiting!!!" << endl;
165 if (fVerbosityLevel >=3)
166 cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
168 fVolumeMediaMap = new TClonesArray("FlukaVolume",1000);
174 if (fVerbosityLevel >=3)
175 cout << "==> TFluka::~TFluka() destructor called." << endl;
177 delete fGeometryManager;
178 fVolumeMediaMap->Delete();
179 delete fVolumeMediaMap;
182 if (fVerbosityLevel >=3)
183 cout << "<== TFluka::~TFluka() destructor called." << endl;
187 //_____________________________________________________________________________
188 // TFluka control methods
189 //____________________________________________________________________________
190 void TFluka::Init() {
191 if (fVerbosityLevel >=3)
192 cout << "==> TFluka::Init() called." << endl;
194 if (fVerbosityLevel >=2)
195 cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
196 << ") in fluka..." << endl;
197 GLOBAL.lfdrtr = true;
199 if (fVerbosityLevel >=2)
200 cout << "\t* Opening file " << fInputFileName << endl;
201 const char* fname = fInputFileName;
202 fluka_openinp(lunin, PASSCHARA(fname));
204 if (fVerbosityLevel >=2)
205 cout << "\t* Calling flukam..." << endl;
208 if (fVerbosityLevel >=2)
209 cout << "\t* Closing file " << fInputFileName << endl;
210 fluka_closeinp(lunin);
212 if (fVerbosityLevel >=3)
213 cout << "<== TFluka::Init() called." << endl;
219 void TFluka::FinishGeometry() {
221 // Build-up table with region to medium correspondance
225 if (fVerbosityLevel >=3)
226 cout << "==> TFluka::FinishGeometry() called." << endl;
228 // fGeometryManager->Ggclos();
230 FGeometryInit* flugg = FGeometryInit::GetInstance();
232 fMediaByRegion = new Int_t[fNVolumes+2];
233 for (Int_t i = 0; i < fNVolumes; i++)
235 FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]);
236 TString volName = vol->GetName();
237 Int_t media = vol->GetMedium();
238 printf("Finish Geometry: volName, media %d %s %d \n", i, volName.Data(), media);
239 strcpy(tmp, volName.Data());
241 flugg->SetMediumFromName(tmp, media);
244 flugg->BuildMediaMap();
246 if (fVerbosityLevel >=3)
247 cout << "<== TFluka::FinishGeometry() called." << endl;
250 void TFluka::BuildPhysics() {
251 if (fVerbosityLevel >=3)
252 cout << "==> TFluka::BuildPhysics() called." << endl;
255 if (fVerbosityLevel >=3)
256 cout << "<== TFluka::BuildPhysics() called." << endl;
259 void TFluka::ProcessEvent() {
260 if (fVerbosityLevel >=3)
261 cout << "==> TFluka::ProcessEvent() called." << endl;
263 if (fVerbosityLevel >=3)
264 cout << "<== TFluka::ProcessEvent() called." << endl;
268 void TFluka::ProcessRun(Int_t nevent) {
269 if (fVerbosityLevel >=3)
270 cout << "==> TFluka::ProcessRun(" << nevent << ") called."
273 if (fVerbosityLevel >=2) {
274 cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
275 cout << "\t* Calling flukam again..." << endl;
277 fApplication->GeneratePrimaries();
278 EPISOR.lsouit = true;
281 if (fVerbosityLevel >=3)
282 cout << "<== TFluka::ProcessRun(" << nevent << ") called."
286 //_____________________________________________________________________________
287 // methods for building/management of geometry
288 //____________________________________________________________________________
289 // functions from GCONS
290 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
291 Float_t &dens, Float_t &radl, Float_t &absl,
292 Float_t* ubuf, Int_t& nbuf) {
294 fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
297 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
298 Double_t &dens, Double_t &radl, Double_t &absl,
299 Double_t* ubuf, Int_t& nbuf) {
301 fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
304 // detector composition
305 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
306 Double_t z, Double_t dens, Double_t radl, Double_t absl,
307 Float_t* buf, Int_t nwbuf) {
310 ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
312 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
313 Double_t z, Double_t dens, Double_t radl, Double_t absl,
314 Double_t* buf, Int_t nwbuf) {
317 ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
320 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
321 Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
324 ->Mixture(kmat, name, a, z, dens, nlmat, wmat);
326 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
327 Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
330 ->Mixture(kmat, name, a, z, dens, nlmat, wmat);
333 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
334 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
335 Double_t stemax, Double_t deemax, Double_t epsil,
336 Double_t stmin, Float_t* ubuf, Int_t nbuf) {
339 ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
340 epsil, stmin, ubuf, nbuf);
342 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
343 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
344 Double_t stemax, Double_t deemax, Double_t epsil,
345 Double_t stmin, Double_t* ubuf, Int_t nbuf) {
348 ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
349 epsil, stmin, ubuf, nbuf);
352 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
353 Double_t thetaY, Double_t phiY, Double_t thetaZ,
357 ->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
360 void TFluka::Gstpar(Int_t itmed, const char *param, Double_t parval) {
362 fGeometryManager->Gstpar(itmed, param, parval);
365 // functions from GGEOM
366 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
367 Float_t *upar, Int_t np) {
369 // fVolumeMediaMap[TString(name)] = nmed;
370 TClonesArray &lvols = *fVolumeMediaMap;
371 new(lvols[fNVolumes++])
372 FlukaVolume(name, nmed);
373 return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
375 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
376 Double_t *upar, Int_t np) {
378 TClonesArray &lvols = *fVolumeMediaMap;
379 new(lvols[fNVolumes++])
380 FlukaVolume(name, nmed);
382 return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
385 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
388 fGeometryManager->Gsdvn(name, mother, ndiv, iaxis);
391 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
392 Int_t iaxis, Double_t c0i, Int_t numed) {
394 TClonesArray &lvols = *fVolumeMediaMap;
395 new(lvols[fNVolumes++])
396 FlukaVolume(name, numed);
397 fGeometryManager->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
400 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
401 Int_t iaxis, Int_t numed, Int_t ndvmx) {
403 TClonesArray &lvols = *fVolumeMediaMap;
404 new(lvols[fNVolumes++])
405 FlukaVolume(name, numed);
406 fGeometryManager->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
409 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
410 Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
412 TClonesArray &lvols = *fVolumeMediaMap;
413 new(lvols[fNVolumes++])
414 FlukaVolume(name, numed);
415 fGeometryManager->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
418 void TFluka::Gsord(const char *name, Int_t iax) {
420 fGeometryManager->Gsord(name, iax);
423 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
424 Double_t x, Double_t y, Double_t z, Int_t irot,
427 fGeometryManager->Gspos(name, nr, mother, x, y, z, irot, konly);
430 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
431 Double_t x, Double_t y, Double_t z, Int_t irot,
432 const char *konly, Float_t *upar, Int_t np) {
434 fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
436 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
437 Double_t x, Double_t y, Double_t z, Int_t irot,
438 const char *konly, Double_t *upar, Int_t np) {
440 fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
443 void TFluka::Gsbool(const char* onlyVolName, const char* manyVolName) {
445 fGeometryManager->Gsbool(onlyVolName, manyVolName);
448 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t *ppckov,
449 Float_t *absco, Float_t *effic, Float_t *rindex) {
451 fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
453 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
454 Double_t *absco, Double_t *effic, Double_t *rindex) {
456 fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
460 void TFluka::WriteEuclid(const char* fileName, const char* topVol,
461 Int_t number, Int_t nlevel) {
463 fGeometryManager->WriteEuclid(fileName, topVol, number, nlevel);
468 //_____________________________________________________________________________
469 // methods needed by the stepping
470 //____________________________________________________________________________
472 Int_t TFluka::GetMedium() const {
473 FGeometryInit* flugg = FGeometryInit::GetInstance();
474 return flugg->GetMedium(fCurrentFlukaRegion);
479 //____________________________________________________________________________
480 // ID <--> PDG transformations
481 //_____________________________________________________________________________
482 Int_t TFluka::IdFromPDG(Int_t pdg) const
485 // Return Fluka code from PDG and pseudo ENDF code
487 // MCIHAD() goes from pdg to fluka internal.
488 Int_t intfluka = mcihad(pdg);
489 // KPTOIP array goes from internal to official
490 return GetFlukaKPTOIP(intfluka);
493 Int_t TFluka::PDGFromId(Int_t id) const
496 // Return PDG code and pseudo ENDF code from Fluka code
498 //IPTOKP array goes from official to internal
499 Int_t intfluka = GetFlukaIPTOKP(id);
500 //MPKDHA() goes from internal to PDG
501 return mpdgha(intfluka);
508 //_____________________________________________________________________________
509 // methods for step management
510 //____________________________________________________________________________
512 // dynamic properties
514 void TFluka::TrackPosition(TLorentzVector& position) const
516 // Return the current position in the master reference frame of the
517 // track being transported
518 // TRACKR.atrack = age of the particle
519 // TRACKR.xtrack = x-position of the last point
520 // TRACKR.ytrack = y-position of the last point
521 // TRACKR.ztrack = z-position of the last point
522 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
523 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
524 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
525 position.SetT(TRACKR.atrack);
528 void TFluka::TrackMomentum(TLorentzVector& momentum) const
530 // Return the direction and the momentum (GeV/c) of the track
531 // currently being transported
532 // TRACKR.ptrack = momentum of the particle (not always defined, if
533 // < 0 must be obtained from etrack)
534 // TRACKR.cx,y,ztrck = direction cosines of the current particle
535 // TRACKR.etrack = total energy of the particle
536 // TRACKR.jtrack = identity number of the particle
537 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
538 if (TRACKR.ptrack >= 0) {
539 momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
540 momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
541 momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
542 momentum.SetE(TRACKR.etrack);
546 Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack]*PAPROP.am[TRACKR.jtrack]);
547 momentum.SetPx(p*TRACKR.cxtrck);
548 momentum.SetPy(p*TRACKR.cytrck);
549 momentum.SetPz(p*TRACKR.cztrck);
550 momentum.SetE(TRACKR.etrack);
555 Double_t TFluka::TrackStep() const
557 // Return the length in centimeters of the current step
558 // TRACKR.ctrack = total curved path
559 return TRACKR.ctrack;
562 Double_t TFluka::TrackLength() const
565 // should be the sum of all steps starting from the beginning of the track
566 // for the time being returns only the length in centimeters of the current step
567 return TRACKR.ctrack;
570 Double_t TFluka::TrackTime() const
572 // Return the current time of flight of the track being transported
573 // TRACKR.atrack = age of the particle
574 return TRACKR.atrack;
577 Double_t TFluka::Edep() const
580 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
581 // -->local energy deposition (the value and the point are not recorded in TRACKR)
582 // but in the variable "rull" of the procedure "endraw.cxx"
583 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
584 // -->no energy loss along the track
585 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
586 // -->energy loss distributed along the track
587 // TRACKR.dtrack = energy deposition of the jth deposition even
588 if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
592 for ( Int_t j=0;j<TRACKR.mtrack;j++) {
593 sum +=TRACKR.dtrack[j];
599 Int_t TFluka::TrackPid() const
601 // Return the id of the particle transported
602 // TRACKR.jtrack = identity number of the particle
603 return PDGFromId(TRACKR.jtrack);
606 Double_t TFluka::TrackCharge() const
608 // Return charge of the track currently transported
609 // PAPROP.ichrge = electric charge of the particle
610 return PAPROP.ichrge[TRACKR.jtrack+6];
613 Double_t TFluka::TrackMass() const
615 // PAPROP.am = particle mass in GeV
616 return PAPROP.am[TRACKR.jtrack+6];
619 Double_t TFluka::Etot() const
621 // TRACKR.etrack = total energy of the particle
622 return TRACKR.etrack;
628 Bool_t TFluka::IsNewTrack() const
631 // True if the track is not at the boundary of the current volume
635 Bool_t TFluka::IsTrackInside() const
637 // True if the track is not at the boundary of the current volume
638 // In Fluka a step is always inside one kind of material
639 // If the step would go behind the region of one material,
640 // it will be shortened to reach only the boundary.
641 // Therefore IsTrackInside() is always true.
645 Bool_t TFluka::IsTrackEntering() const
647 // True if this is the first step of the track in the current volume
648 // Boundary- (X) crossing
649 // Icode = 19: boundary crossing - call from Kaskad
650 // Icode = 29: boundary crossing - call from Emfsco
651 // Icode = 39: boundary crossing - call from Kasneu
652 // Icode = 49: boundary crossing - call from Kashea
653 // Icode = 59: boundary crossing - call from Kasoph
658 fIcode == 59) return 1;
662 Bool_t TFluka::IsTrackExiting() const
664 // True if this is the last step of the track in the current volume
665 // Boundary- (X) crossing
666 // Icode = 19: boundary crossing - call from Kaskad
667 // Icode = 29: boundary crossing - call from Emfsco
668 // Icode = 39: boundary crossing - call from Kasneu
669 // Icode = 49: boundary crossing - call from Kashea
670 // Icode = 59: boundary crossing - call from Kasoph
675 fIcode == 59) return 1;
679 Bool_t TFluka::IsTrackOut() const
681 // True if the track is out of the setup
683 // Icode = 14: escape - call from Kaskad
684 // Icode = 23: escape - call from Emfsco
685 // Icode = 32: escape - call from Kasneu
686 // Icode = 40: escape - call from Kashea
687 // Icode = 51: escape - call from Kasoph
692 fIcode == 51) return 1;
696 Bool_t TFluka::IsTrackDisappeared() const
698 // means all inelastic interactions and decays
699 // Icode = 11: inelastic interaction recoil - call from Kaskad
700 if (fIcode == 11) return 1;
704 Bool_t TFluka::IsTrackStop() const
706 // True if the track energy has fallen below the threshold
707 // means stopped by signal or below energy threshold
708 // Icode = 12: stopping particle - call from Kaskad
709 // Icode = 15: time kill - call from Kaskad
710 // Icode = 21: below threshold, iarg=1 - call from Emfsco
711 // Icode = 22: below threshold, iarg=2 - call from Emfsco
712 // Icode = 24: time kill - call from Emfsco
713 // Icode = 31: below threshold - call from Kasneu
714 // Icode = 33: time kill - call from Kasneu
715 // Icode = 41: time kill - call from Kashea
716 // Icode = 52: time kill - call from Kasoph
725 fIcode == 52) return 1;
729 Bool_t TFluka::IsTrackAlive() const
731 // means not disappeared or not out
732 if (IsTrackDisappeared() || IsTrackOut() ) return 0;
740 //Int_t NSecondaries() const
741 // Number of secondary particles generated in the current step
746 //void GetSecondary(Int_t isec, Int_t& particleId, TLorentzVector& position,
747 // TLorentzVector& momentum)
750 // will come from FINUC when called from USDRAW
753 //TMCProcess ProdProcess(Int_t isec) const
754 // Name of the process that has produced the secondary particles
755 // in the current step
757 // will come from FINUC when called from USDRAW
760 //Int_t StepProcesses(TArrayI &proc) const
761 // Return processes active in the current step
763 //ck = total energy of the particl ????????????????
767 void TFluka::FutoTest() const
769 Int_t icode, mreg, newreg, medium;
770 Double_t rull, xsco, ysco, zsco;
775 cout << " icode=" << icode << endl;
776 else if (icode > 0 && icode <= 5) {
778 medium = GetMedium();
779 cout << " icode=" << icode
781 << " medium=" << medium
784 else if (icode >= 10 && icode <= 15) {
786 medium = GetMedium();
791 cout << " icode=" << icode
793 << " medium=" << medium
797 << " zsco=" << zsco << endl;
799 else if (icode >= 20 && icode <= 24) {
801 medium = GetMedium();
806 cout << " icode=" << icode
808 << " medium=" << medium
812 << " zsco=" << zsco << endl;
814 else if (icode >= 30 && icode <= 33) {
816 medium = GetMedium();
821 cout << " icode=" << icode
823 << " medium=" << medium
827 << " zsco=" << zsco << endl;
829 else if (icode >= 40 && icode <= 41) {
831 medium = GetMedium();
836 cout << " icode=" << icode
838 << " medium=" << medium
842 << " zsco=" << zsco << endl;
844 else if (icode >= 50 && icode <= 52) {
846 medium = GetMedium();
851 cout << " icode=" << icode
853 << " medium=" << medium
857 << " zsco=" << zsco << endl;
859 else if (icode >= 100 && icode <= 105) {
861 medium = GetMedium();
865 cout << " icode=" << icode
867 << " medium=" << medium
870 << " zsco=" << zsco << endl;
872 else if (icode == 208) {
874 medium = GetMedium();
878 cout << " icode=" << icode
880 << " medium=" << medium
883 << " zsco=" << zsco << endl;
885 else if (icode == 210) {
887 medium = GetMedium();
891 cout << " icode=" << icode
893 << " medium=" << medium
896 << " zsco=" << zsco << endl;
898 else if (icode == 212) {
900 medium = GetMedium();
904 cout << " icode=" << icode
906 << " medium=" << medium
909 << " zsco=" << zsco << endl;
911 else if (icode >= 214 && icode <= 215) {
913 medium = GetMedium();
917 cout << " icode=" << icode
919 << " medium=" << medium
922 << " zsco=" << zsco << endl;
924 else if (icode == 217) {
926 medium = GetMedium();
930 cout << " icode=" << icode
932 << " medium=" << medium
935 << " zsco=" << zsco << endl;
937 else if (icode == 219) {
939 medium = GetMedium();
943 cout << " icode=" << icode
945 << " medium=" << medium
948 << " zsco=" << zsco << endl;
950 else if (icode == 221) {
952 medium = GetMedium();
956 cout << " icode=" << icode
958 << " medium=" << medium
961 << " zsco=" << zsco << endl;
963 else if (icode == 225) {
965 medium = GetMedium();
969 cout << " icode=" << icode
971 << " medium=" << medium
974 << " zsco=" << zsco << endl;
976 else if (icode == 300) {
978 medium = GetMedium();
982 cout << " icode=" << icode
984 << " medium=" << medium
987 << " zsco=" << zsco << endl;
989 else if (icode == 400) {
991 medium = GetMedium();
995 cout << " icode=" << icode
997 << " medium=" << medium
1000 << " zsco=" << zsco << endl;
1002 else if (icode == 19) {
1004 medium = GetMedium();
1005 newreg = GetNewreg();
1009 cout << " icode=" << icode
1011 << " medium=" << medium
1012 << " newreg=" << newreg
1015 << " zsco=" << zsco << endl;
1017 else if (icode == 29) {
1019 medium = GetMedium();
1020 newreg = GetNewreg();
1024 cout << " icode=" << icode
1026 << " medium=" << medium
1027 << " newreg=" << newreg
1030 << " zsco=" << zsco << endl;
1032 else if (icode == 39) {
1034 medium = GetMedium();
1035 newreg = GetNewreg();
1039 cout << " icode=" << icode
1041 << " medium=" << medium
1042 << " newreg=" << newreg
1045 << " zsco=" << zsco << endl;
1047 else if (icode == 49) {
1049 medium = GetMedium();
1050 newreg = GetNewreg();
1054 cout << " icode=" << icode
1056 << " medium=" << medium
1057 << " newreg=" << newreg
1060 << " zsco=" << zsco << endl;
1062 else if (icode == 59) {
1064 medium = GetMedium();
1065 newreg = GetNewreg();
1069 cout << " icode=" << icode
1071 << " medium=" << medium
1072 << " newreg=" << newreg
1075 << " zsco=" << zsco << endl;
1079 } // end of FutoTest
1085 //_____________________________________________________________________________
1086 // methods for step management
1087 //____________________________________________________________________________
1089 // dynamic properties
1091 void TFluka::TrackPosition(TLorentzVector& position) const
1093 // Return the current position in the master reference frame of the
1094 // track being transported
1095 // TRACKR.atrack = age of the particle
1096 // TRACKR.xtrack = x-position of the last point
1097 // TRACKR.ytrack = y-position of the last point
1098 // TRACKR.ztrack = z-position of the last point
1099 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1100 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1101 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1102 position.SetT(TRACKR.atrack);
1105 void TFluka::TrackMomentum(TLorentzVector& momentum) const
1107 // Return the direction and the momentum (GeV/c) of the track
1108 // currently being transported
1109 // TRACKR.ptrack = momentum of the particle (not always defined, if
1110 // < 0 must be obtained from etrack)
1111 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1112 // TRACKR.etrack = total energy of the particle
1113 // TRACKR.jtrack = identity number of the particle
1114 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1115 if (TRACKR.ptrack >= 0) {
1116 momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
1117 momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
1118 momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
1119 momentum.SetE(TRACKR.etrack);
1123 Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack]*PAPROP.am[TRACKR.jtrack]);
1124 momentum.SetPx(p*TRACKR.cxtrck);
1125 momentum.SetPy(p*TRACKR.cytrck);
1126 momentum.SetPz(p*TRACKR.cztrck);
1127 momentum.SetE(TRACKR.etrack);
1132 Double_t TFluka::TrackStep() const
1134 // Return the length in centimeters of the current step
1135 // TRACKR.ctrack = total curved path
1136 return TRACKR.ctrack;
1139 Double_t TFluka::TrackLength() const
1142 // should be the sum of all steps starting from the beginning of the track
1143 // for the time being returns only the length in centimeters of the current step
1144 return TRACKR.ctrack;
1147 Double_t TFluka::TrackTime() const
1149 // Return the current time of flight of the track being transported
1150 // TRACKR.atrack = age of the particle
1151 return TRACKR.atrack;
1154 Double_t TFluka::Edep() const
1156 // Energy deposition
1157 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
1158 // -->local energy deposition (the value and the point are not recorded in TRACKR)
1159 // but in the variable "rull" of the procedure "endraw.cxx"
1160 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
1161 // -->no energy loss along the track
1162 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
1163 // -->energy loss distributed along the track
1164 // TRACKR.dtrack = energy deposition of the jth deposition even
1165 if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
1169 for ( Int_t j=0;j<TRACKR.mtrack;j++) {
1170 sum +=TRACKR.dtrack[j];
1176 Int_t TFluka::TrackPid() const
1178 // Return the id of the particle transported
1179 // TRACKR.jtrack = identity number of the particle
1180 return PDGFromId(TRACKR.jtrack);
1183 Double_t TFluka::TrackCharge() const
1185 // Return charge of the track currently transported
1186 // PAPROP.ichrge = electric charge of the particle
1187 return PAPROP.ichrge[TRACKR.jtrack+6];
1190 Double_t TFluka::TrackMass() const
1192 // PAPROP.am = particle mass in GeV
1193 return PAPROP.am[TRACKR.jtrack+6];
1196 Double_t TFluka::Etot() const
1198 // TRACKR.etrack = total energy of the particle
1199 return TRACKR.etrack;
1205 Bool_t TFluka::IsNewTrack() const
1208 // True if the track is not at the boundary of the current volume
1209 // Not true in some cases in bxdraw - to be solved
1213 Bool_t TFluka::IsTrackInside() const
1215 // True if the track is not at the boundary of the current volume
1216 // In Fluka a step is always inside one kind of material
1217 // If the step would go behind the region of one material,
1218 // it will be shortened to reach only the boundary.
1219 // Therefore IsTrackInside() is always true.
1220 // Not true in some cases in bxdraw - to be solved
1224 Bool_t TFluka::IsTrackEntering() const
1226 // True if this is the first step of the track in the current volume
1227 // Boundary- (X) crossing
1228 // Icode = 19: boundary crossing - call from Kaskad
1229 // Icode = 29: boundary crossing - call from Emfsco
1230 // Icode = 39: boundary crossing - call from Kasneu
1231 // Icode = 49: boundary crossing - call from Kashea
1232 // Icode = 59: boundary crossing - call from Kasoph
1237 fIcode == 59) return 1;
1241 Bool_t TFluka::IsTrackExiting() const
1243 // True if this is the last step of the track in the current volume
1244 // Boundary- (X) crossing
1245 // Icode = 19: boundary crossing - call from Kaskad
1246 // Icode = 29: boundary crossing - call from Emfsco
1247 // Icode = 39: boundary crossing - call from Kasneu
1248 // Icode = 49: boundary crossing - call from Kashea
1249 // Icode = 59: boundary crossing - call from Kasoph
1254 fIcode == 59) return 1;
1258 Bool_t TFluka::IsTrackOut() const
1260 // True if the track is out of the setup
1262 // Icode = 14: escape - call from Kaskad
1263 // Icode = 23: escape - call from Emfsco
1264 // Icode = 32: escape - call from Kasneu
1265 // Icode = 40: escape - call from Kashea
1266 // Icode = 51: escape - call from Kasoph
1271 fIcode == 51) return 1;
1275 Bool_t TFluka::IsTrackDisappeared() const
1277 // means all inelastic interactions and decays
1278 // fIcode from usdraw
1279 if (fIcode == 101 || // inelastic interaction
1280 fIcode == 102 || // particle decay
1281 fIcode == 214 || // in-flight annihilation
1282 fIcode == 215 || // annihilation at rest
1283 fIcode == 217 || // pair production
1284 fIcode == 221) return 1;
1288 Bool_t TFluka::IsTrackStop() const
1290 // True if the track energy has fallen below the threshold
1291 // means stopped by signal or below energy threshold
1292 // Icode = 12: stopping particle - call from Kaskad
1293 // Icode = 15: time kill - call from Kaskad
1294 // Icode = 21: below threshold, iarg=1 - call from Emfsco
1295 // Icode = 22: below threshold, iarg=2 - call from Emfsco
1296 // Icode = 24: time kill - call from Emfsco
1297 // Icode = 31: below threshold - call from Kasneu
1298 // Icode = 33: time kill - call from Kasneu
1299 // Icode = 41: time kill - call from Kashea
1300 // Icode = 52: time kill - call from Kasoph
1309 fIcode == 52) return 1;
1313 Bool_t TFluka::IsTrackAlive() const
1315 // means not disappeared or not out
1316 if (IsTrackDisappeared() || IsTrackOut() ) return 0;
1324 Int_t TFluka::NSecondaries() const
1325 // Number of secondary particles generated in the current step
1326 // fIcode = 100 = elastic interaction
1327 // fIcode = 101 = inelastic interaction
1328 // fIcode = 102 = particle decay
1329 // fIcode = 103 = delta ray
1330 // fIcode = 104 = pair production
1331 // fIcode = 105 = bremsstrahlung
1333 if (fIcode >= 100 && fIcode <= 105)
1334 return FINUC.np + FHEAVY.npheav;
1339 void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
1340 TLorentzVector& position, TLorentzVector& momentum)
1342 // fIcode = 100 = elastic interaction
1343 // fIcode = 101 = inelastic interaction
1344 // fIcode = 102 = particle decay
1345 // fIcode = 103 = delta ray
1346 // fIcode = 104 = pair production
1347 // fIcode = 105 = bremsstrahlung
1349 if (fIcode >= 100 && fIcode <= 105) {
1350 if (isec >= 0 && isec < FINUC.np) {
1351 particleId = PDGFromId(FINUC.kpart[isec]);
1352 position.SetX(fXsco);
1353 position.SetY(fYsco);
1354 position.SetZ(fZsco);
1355 position.SetT(FINUC.agesec[isec]);
1356 momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
1357 momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
1358 momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
1359 momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
1361 if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
1362 Int_t jsec = isec - FINUC.np;
1363 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
1364 position.SetX(fXsco);
1365 position.SetY(fYsco);
1366 position.SetZ(fZsco);
1367 position.SetT(FHEAVY.agheav[jsec]);
1368 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
1369 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
1370 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
1371 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
1372 momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
1373 else if (FHEAVY.tkheav[jsec] > 6)
1374 momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
1379 //TMCProcess ProdProcess(Int_t isec) const
1380 // Name of the process that has produced the secondary particles
1381 // in the current step
1383 // will come from FINUC when called from USDRAW
1386 //Int_t StepProcesses(TArrayI &proc) const
1387 // Return processes active in the current step
1389 //ck = total energy of the particl ????????????????
1393 // ===============================================================
1394 void TFluka::FutoTest()
1396 Int_t icode, mreg, newreg, particleId;
1398 Double_t rull, xsco, ysco, zsco;
1399 TLorentzVector position, momentum;
1402 cout << " icode=" << icode << endl;
1404 cout << "TLorentzVector positionX=" << position.X()
1405 << "positionY=" << position.Y()
1406 << "positionZ=" << position.Z()
1407 << "timeT=" << position.T() << endl;
1408 cout << "TLorentzVector momentumX=" << momentum.X()
1409 << "momentumY=" << momentum.Y()
1410 << "momentumZ=" << momentum.Z()
1411 << "energyE=" << momentum.E() << endl;
1412 cout << "TrackPid=" << TrackPid() << endl;
1416 else if (icode > 0 && icode <= 5) {
1419 // medium = GetMedium();
1420 cout << " icode=" << icode
1422 // << " medium=" << medium
1424 TrackPosition(position);
1425 TrackMomentum(momentum);
1426 cout << "TLorentzVector positionX=" << position.X()
1427 << "positionY=" << position.Y()
1428 << "positionZ=" << position.Z()
1429 << "timeT=" << position.T() << endl;
1430 cout << "TLorentzVector momentumX=" << momentum.X()
1431 << "momentumY=" << momentum.Y()
1432 << "momentumZ=" << momentum.Z()
1433 << "energyE=" << momentum.E() << endl;
1434 cout << "TrackStep=" << TrackStep() << endl;
1435 cout << "TrackLength=" << TrackLength() << endl;
1436 cout << "TrackTime=" << TrackTime() << endl;
1437 cout << "Edep=" << Edep() << endl;
1438 cout << "TrackPid=" << TrackPid() << endl;
1439 cout << "TrackCharge=" << TrackCharge() << endl;
1440 cout << "TrackMass=" << TrackMass() << endl;
1441 cout << "Etot=" << Etot() << endl;
1442 cout << "IsNewTrack=" << IsNewTrack() << endl;
1443 cout << "IsTrackInside=" << IsTrackInside() << endl;
1444 cout << "IsTrackEntering=" << IsTrackEntering() << endl;
1445 cout << "IsTrackExiting=" << IsTrackExiting() << endl;
1446 cout << "IsTrackOut=" << IsTrackOut() << endl;
1447 cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
1448 cout << "IsTrackAlive=" << IsTrackAlive() << endl;
1451 else if((icode >= 10 && icode <= 15) ||
1452 (icode >= 20 && icode <= 24) ||
1453 (icode >= 30 && icode <= 33) ||
1454 (icode >= 40 && icode <= 41) ||
1455 (icode >= 50 && icode <= 52)) {
1458 // medium = GetMedium();
1463 cout << " icode=" << icode
1465 // << " medium=" << medium
1469 << " zsco=" << zsco << endl;
1470 TrackPosition(position);
1471 TrackMomentum(momentum);
1472 cout << "Edep=" << Edep() << endl;
1473 cout << "Etot=" << Etot() << endl;
1474 cout << "TrackPid=" << TrackPid() << endl;
1475 cout << "TrackCharge=" << TrackCharge() << endl;
1476 cout << "TrackMass=" << TrackMass() << endl;
1477 cout << "IsTrackOut=" << IsTrackOut() << endl;
1478 cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
1479 cout << "IsTrackStop=" << IsTrackStop() << endl;
1480 cout << "IsTrackAlive=" << IsTrackAlive() << endl;
1483 else if((icode >= 100 && icode <= 105) ||
1487 (icode >= 214 && icode <= 215) ||
1496 // medium = GetMedium();
1500 cout << " icode=" << icode
1502 // << " medium=" << medium
1505 << " zsco=" << zsco << endl;
1506 cout << "TrackPid=" << TrackPid() << endl;
1507 cout << "NSecondaries=" << NSecondaries() << endl;
1508 for (Int_t isec=0; isec< NSecondaries(); isec++) {
1509 //void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
1510 // TLorentzVector& position, TLorentzVector& momentum)
1511 TFluka::GetSecondary(isec, particleId, position, momentum);
1512 cout << "TLorentzVector positionX=" << position.X()
1513 << "positionY=" << position.Y()
1514 << "positionZ=" << position.Z()
1515 << "timeT=" << position.T() << endl;
1516 cout << "TLorentzVector momentumX=" << momentum.X()
1517 << "momentumY=" << momentum.Y()
1518 << "momentumZ=" << momentum.Z()
1519 << "energyE=" << momentum.E() << endl;
1520 cout << "TrackPid=" << particleId << endl;
1525 else if((icode == 19) ||
1531 // medium = GetMedium();
1532 newreg = GetNewreg();
1536 cout << " icode=" << icode
1538 // << " medium=" << medium
1539 << " newreg=" << newreg
1542 << " zsco=" << zsco << endl;
1545 // ====================================================================
1550 } // end of FutoTest