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.12 2003/01/31 12:18:53 morsch
19 Corrected indices. (E. Futo)
21 Revision 1.9 2002/12/06 12:41:29 morsch
22 Mess from last merge cleaned up.
24 Revision 1.8 2002/12/06 12:28:44 morsch
25 Region to media mapping corrected and improved.
27 Revision 1.7 2002/12/06 12:21:32 morsch
28 User stepping methods added (E. Futo)
30 Revision 1.6 2002/11/21 18:40:06 iglez2
33 Revision 1.5 2002/11/07 17:59:10 iglez2
34 Included the geometry through geant4_vmc/FLUGG
36 Revision 1.4 2002/11/04 16:00:46 iglez2
37 The conversion between ID and PDG now uses Fluka routines and arrays which is more consistent.
39 Revision 1.3 2002/10/22 15:12:14 alibrary
40 Introducing Riostream.h
42 Revision 1.2 2002/10/14 14:57:40 hristov
43 Merging the VirtualMC branch to the main development branch (HEAD)
45 Revision 1.1.2.8 2002/10/08 16:33:17 iglez2
46 LSOUIT is set to true before the second call to flukam.
48 Revision 1.1.2.7 2002/10/08 09:30:37 iglez2
49 Solved stupid missing ;
51 Revision 1.1.2.6 2002/10/07 13:40:22 iglez2
52 First implementations of the PDG <--> Fluka Id conversion routines
54 Revision 1.1.2.5 2002/09/26 16:26:03 iglez2
56 Call to gAlice->Generator()->Generate()
58 Revision 1.1.2.4 2002/09/26 13:22:23 iglez2
59 Naive implementation of ProcessRun and ProcessEvent
60 Opening/Closing of input file (fInputFileName) with FORTRAN unit 5 before/after the first call to flukam inside Init()
62 Revision 1.1.2.3 2002/09/20 15:35:51 iglez2
63 Modification of LFDRTR. Value is passed to FLUKA !!!
65 Revision 1.1.2.2 2002/09/18 14:34:44 iglez2
66 Revised version with all pure virtual methods implemented
68 Revision 1.1.2.1 2002/07/24 08:49:41 alibrary
69 Adding TFluka to VirtualMC
71 Revision 1.1 2002/07/05 13:10:07 morsch
72 First commit of Fluka interface.
76 #include <Riostream.h>
78 #include "TClonesArray.h"
80 #include "TCallf77.h" //For the fortran calls
81 #include "Fdblprc.h" //(DBLPRC) fluka common
82 #include "Fepisor.h" //(EPISOR) fluka common
83 #include "Ffinuc.h" //(FINUC) fluka common
84 #include "Fiounit.h" //(IOUNIT) fluka common
85 #include "Fpaprop.h" //(PAPROP) fluka common
86 #include "Fpart.h" //(PART) fluka common
87 #include "Ftrackr.h" //(TRACKR) fluka common
88 #include "Fpaprop.h" //(PAPROP) fluka common
89 #include "Ffheavy.h" //(FHEAVY) fluka common
91 #include "TVirtualMC.h"
92 #include "TG4GeometryManager.h" //For the geometry management
93 #include "TG4DetConstruction.h" //For the detector construction
95 #include "FGeometryInit.hh"
96 #include "TLorentzVector.h"
97 #include "FlukaVolume.h"
99 // Fluka methods that may be needed.
101 # define flukam flukam_
102 # define fluka_openinp fluka_openinp_
103 # define fluka_closeinp fluka_closeinp_
104 # define mcihad mcihad_
105 # define mpdgha mpdgha_
107 # define flukam FLUKAM
108 # define fluka_openinp FLUKA_OPENINP
109 # define fluka_closeinp FLUKA_CLOSEINP
110 # define mcihad MCIHAD
111 # define mpdgha MPDGHA
117 // Prototypes for FLUKA functions
119 void type_of_call flukam(const int&);
120 void type_of_call fluka_openinp(const int&, DEFCHARA);
121 void type_of_call fluka_closeinp(const int&);
122 int type_of_call mcihad(const int&);
123 int type_of_call mpdgha(const int&);
127 // Class implementation for ROOT
132 //----------------------------------------------------------------------------
133 // TFluka constructors and destructors.
134 //____________________________________________________________________________
140 fCurrentFlukaRegion(-1)
143 // Default constructor
147 TFluka::TFluka(const char *title, Int_t verbosity)
148 :TVirtualMC("TFluka",title),
149 fVerbosityLevel(verbosity),
152 fCurrentFlukaRegion(-1)
154 if (fVerbosityLevel >=3)
155 cout << "==> TFluka::TFluka(" << title << ") constructor called." << endl;
158 // create geometry manager
159 if (fVerbosityLevel >=2)
160 cout << "\t* Creating G4 Geometry manager..." << endl;
161 fGeometryManager = new TG4GeometryManager();
162 if (fVerbosityLevel >=2)
163 cout << "\t* Creating G4 Detector..." << endl;
164 fDetector = new TG4DetConstruction();
165 FGeometryInit* geominit = FGeometryInit::GetInstance();
167 geominit->setDetConstruction(fDetector);
169 cerr << "ERROR: Could not create FGeometryInit!" << endl;
170 cerr << " Exiting!!!" << endl;
174 if (fVerbosityLevel >=3)
175 cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
177 fVolumeMediaMap = new TClonesArray("FlukaVolume",1000);
183 if (fVerbosityLevel >=3)
184 cout << "==> TFluka::~TFluka() destructor called." << endl;
186 delete fGeometryManager;
187 fVolumeMediaMap->Delete();
188 delete fVolumeMediaMap;
191 if (fVerbosityLevel >=3)
192 cout << "<== TFluka::~TFluka() destructor called." << endl;
196 //_____________________________________________________________________________
197 // TFluka control methods
198 //____________________________________________________________________________
199 void TFluka::Init() {
200 if (fVerbosityLevel >=3)
201 cout << "==> TFluka::Init() called." << endl;
203 if (fVerbosityLevel >=2)
204 cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
205 << ") in fluka..." << endl;
206 GLOBAL.lfdrtr = true;
208 if (fVerbosityLevel >=2)
209 cout << "\t* Opening file " << fInputFileName << endl;
210 const char* fname = fInputFileName;
211 fluka_openinp(lunin, PASSCHARA(fname));
213 if (fVerbosityLevel >=2)
214 cout << "\t* Calling flukam..." << endl;
217 if (fVerbosityLevel >=2)
218 cout << "\t* Closing file " << fInputFileName << endl;
219 fluka_closeinp(lunin);
221 if (fVerbosityLevel >=3)
222 cout << "<== TFluka::Init() called." << endl;
228 void TFluka::FinishGeometry() {
230 // Build-up table with region to medium correspondance
234 if (fVerbosityLevel >=3)
235 cout << "==> TFluka::FinishGeometry() called." << endl;
237 // fGeometryManager->Ggclos();
239 FGeometryInit* flugg = FGeometryInit::GetInstance();
241 fMediaByRegion = new Int_t[fNVolumes+2];
242 for (Int_t i = 0; i < fNVolumes; i++)
244 FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]);
245 TString volName = vol->GetName();
246 Int_t media = vol->GetMedium();
247 printf("Finish Geometry: volName, media %d %s %d \n", i, volName.Data(), media);
248 strcpy(tmp, volName.Data());
250 flugg->SetMediumFromName(tmp, media, i+1);
251 fMediaByRegion[i] = media;
254 flugg->BuildMediaMap();
256 if (fVerbosityLevel >=3)
257 cout << "<== TFluka::FinishGeometry() called." << endl;
260 void TFluka::BuildPhysics() {
261 if (fVerbosityLevel >=3)
262 cout << "==> TFluka::BuildPhysics() called." << endl;
265 if (fVerbosityLevel >=3)
266 cout << "<== TFluka::BuildPhysics() called." << endl;
269 void TFluka::ProcessEvent() {
270 if (fVerbosityLevel >=3)
271 cout << "==> TFluka::ProcessEvent() called." << endl;
272 fApplication->GeneratePrimaries();
273 EPISOR.lsouit = true;
275 if (fVerbosityLevel >=3)
276 cout << "<== TFluka::ProcessEvent() called." << endl;
280 void TFluka::ProcessRun(Int_t nevent) {
281 if (fVerbosityLevel >=3)
282 cout << "==> TFluka::ProcessRun(" << nevent << ") called."
285 if (fVerbosityLevel >=2) {
286 cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
287 cout << "\t* Calling flukam again..." << endl;
289 fApplication->InitGeometry();
290 fApplication->BeginEvent();
292 fApplication->FinishEvent();
293 if (fVerbosityLevel >=3)
294 cout << "<== TFluka::ProcessRun(" << nevent << ") called."
299 //_____________________________________________________________________________
300 // methods for building/management of geometry
301 //____________________________________________________________________________
302 // functions from GCONS
303 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
304 Float_t &dens, Float_t &radl, Float_t &absl,
305 Float_t* ubuf, Int_t& nbuf) {
307 fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
310 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
311 Double_t &dens, Double_t &radl, Double_t &absl,
312 Double_t* ubuf, Int_t& nbuf) {
314 fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
317 // detector composition
318 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
319 Double_t z, Double_t dens, Double_t radl, Double_t absl,
320 Float_t* buf, Int_t nwbuf) {
323 ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
325 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
326 Double_t z, Double_t dens, Double_t radl, Double_t absl,
327 Double_t* buf, Int_t nwbuf) {
330 ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
333 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
334 Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
337 ->Mixture(kmat, name, a, z, dens, nlmat, wmat);
339 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
340 Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
343 ->Mixture(kmat, name, a, z, dens, nlmat, wmat);
346 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
347 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
348 Double_t stemax, Double_t deemax, Double_t epsil,
349 Double_t stmin, Float_t* ubuf, Int_t nbuf) {
352 ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
353 epsil, stmin, ubuf, nbuf);
355 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
356 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
357 Double_t stemax, Double_t deemax, Double_t epsil,
358 Double_t stmin, Double_t* ubuf, Int_t nbuf) {
361 ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
362 epsil, stmin, ubuf, nbuf);
365 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
366 Double_t thetaY, Double_t phiY, Double_t thetaZ,
370 ->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
373 void TFluka::Gstpar(Int_t itmed, const char *param, Double_t parval) {
375 fGeometryManager->Gstpar(itmed, param, parval);
378 // functions from GGEOM
379 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
380 Float_t *upar, Int_t np) {
382 // fVolumeMediaMap[TString(name)] = nmed;
383 printf("TFluka::Gsvolu() name = %s, nmed = %d\n", name, nmed);
385 TClonesArray &lvols = *fVolumeMediaMap;
386 new(lvols[fNVolumes++])
387 FlukaVolume(name, nmed);
388 return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
390 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
391 Double_t *upar, Int_t np) {
393 TClonesArray &lvols = *fVolumeMediaMap;
394 new(lvols[fNVolumes++])
395 FlukaVolume(name, nmed);
397 return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
400 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
403 // The medium of the daughter is the one of the mother
404 Int_t volid = TFluka::VolId(mother);
405 Int_t med = TFluka::VolId2Mate(volid);
406 TClonesArray &lvols = *fVolumeMediaMap;
407 new(lvols[fNVolumes++])
408 FlukaVolume(name, med);
409 fGeometryManager->Gsdvn(name, mother, ndiv, iaxis);
412 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
413 Int_t iaxis, Double_t c0i, Int_t numed) {
415 TClonesArray &lvols = *fVolumeMediaMap;
416 new(lvols[fNVolumes++])
417 FlukaVolume(name, numed);
418 fGeometryManager->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
421 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
422 Int_t iaxis, Int_t numed, Int_t ndvmx) {
424 TClonesArray &lvols = *fVolumeMediaMap;
425 new(lvols[fNVolumes++])
426 FlukaVolume(name, numed);
427 fGeometryManager->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
430 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
431 Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
433 TClonesArray &lvols = *fVolumeMediaMap;
434 new(lvols[fNVolumes++])
435 FlukaVolume(name, numed);
436 fGeometryManager->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
439 void TFluka::Gsord(const char *name, Int_t iax) {
441 fGeometryManager->Gsord(name, iax);
444 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
445 Double_t x, Double_t y, Double_t z, Int_t irot,
448 fGeometryManager->Gspos(name, nr, mother, x, y, z, irot, konly);
451 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
452 Double_t x, Double_t y, Double_t z, Int_t irot,
453 const char *konly, Float_t *upar, Int_t np) {
455 fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
457 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
458 Double_t x, Double_t y, Double_t z, Int_t irot,
459 const char *konly, Double_t *upar, Int_t np) {
461 fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
464 void TFluka::Gsbool(const char* onlyVolName, const char* manyVolName) {
466 fGeometryManager->Gsbool(onlyVolName, manyVolName);
469 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t *ppckov,
470 Float_t *absco, Float_t *effic, Float_t *rindex) {
472 fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
474 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
475 Double_t *absco, Double_t *effic, Double_t *rindex) {
477 fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
481 void TFluka::WriteEuclid(const char* fileName, const char* topVol,
482 Int_t number, Int_t nlevel) {
484 fGeometryManager->WriteEuclid(fileName, topVol, number, nlevel);
489 //_____________________________________________________________________________
490 // methods needed by the stepping
491 //____________________________________________________________________________
493 Int_t TFluka::GetMedium() const {
495 // Get the medium number for the current fluka region
497 FGeometryInit* flugg = FGeometryInit::GetInstance();
498 return flugg->GetMedium(fCurrentFlukaRegion);
503 //____________________________________________________________________________
504 // ID <--> PDG transformations
505 //_____________________________________________________________________________
506 Int_t TFluka::IdFromPDG(Int_t pdg) const
509 // Return Fluka code from PDG and pseudo ENDF code
511 // MCIHAD() goes from pdg to fluka internal.
512 Int_t intfluka = mcihad(pdg);
513 // KPTOIP array goes from internal to official
514 return GetFlukaKPTOIP(intfluka);
517 Int_t TFluka::PDGFromId(Int_t id) const
520 // Return PDG code and pseudo ENDF code from Fluka code
522 //IPTOKP array goes from official to internal
524 printf("PDGFromId: Error id = 0");
528 Int_t intfluka = GetFlukaIPTOKP(id);
530 printf("PDGFromId: Error intfluka = 0");
534 //MPKDHA() goes from internal to PDG
535 return mpdgha(intfluka);
538 //_____________________________________________________________________________
539 // methods for step management
540 //____________________________________________________________________________
544 void TFluka::SetMaxStep(Double_t)
546 // SetMaxStep is dummy procedure in TFluka !
547 cout << "SetMaxStep is dummy procedure in TFluka !" << endl;
550 void TFluka::SetMaxNStep(Int_t)
552 // SetMaxNStep is dummy procedure in TFluka !
553 cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
556 void TFluka::SetUserDecay(Int_t)
558 // SetUserDecay is dummy procedure in TFluka !
559 cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
563 // dynamic properties
565 void TFluka::TrackPosition(TLorentzVector& position) const
567 // Return the current position in the master reference frame of the
568 // track being transported
569 // TRACKR.atrack = age of the particle
570 // TRACKR.xtrack = x-position of the last point
571 // TRACKR.ytrack = y-position of the last point
572 // TRACKR.ztrack = z-position of the last point
573 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
574 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
575 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
576 position.SetT(TRACKR.atrack);
579 void TFluka::TrackMomentum(TLorentzVector& momentum) const
581 // Return the direction and the momentum (GeV/c) of the track
582 // currently being transported
583 // TRACKR.ptrack = momentum of the particle (not always defined, if
584 // < 0 must be obtained from etrack)
585 // TRACKR.cx,y,ztrck = direction cosines of the current particle
586 // TRACKR.etrack = total energy of the particle
587 // TRACKR.jtrack = identity number of the particle
588 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
589 if (TRACKR.ptrack >= 0) {
590 momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
591 momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
592 momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
593 momentum.SetE(TRACKR.etrack);
597 Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
598 momentum.SetPx(p*TRACKR.cxtrck);
599 momentum.SetPy(p*TRACKR.cytrck);
600 momentum.SetPz(p*TRACKR.cztrck);
601 momentum.SetE(TRACKR.etrack);
606 Double_t TFluka::TrackStep() const
608 // Return the length in centimeters of the current step
609 // TRACKR.ctrack = total curved path
610 return TRACKR.ctrack;
613 Double_t TFluka::TrackLength() const
616 // This is the sum of substeps !!!
617 // TRACKR.ctrack = total curved path of the current step
618 // Sum of the substeps is identical to TRACKR.ctrack if the is no mag. field
619 // The sum of all step length starting from the beginning of the track
620 // for the time being returns only the length in centimeters of the current step
622 for ( Int_t j=0;j<TRACKR.ntrack;j++) {
623 sum +=TRACKR.ttrack[j];
628 Double_t TFluka::TrackTime() const
630 // Return the current time of flight of the track being transported
631 // TRACKR.atrack = age of the particle
632 return TRACKR.atrack;
635 Double_t TFluka::Edep() const
638 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
639 // -->local energy deposition (the value and the point are not recorded in TRACKR)
640 // but in the variable "rull" of the procedure "endraw.cxx"
641 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
642 // -->no energy loss along the track
643 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
644 // -->energy loss distributed along the track
645 // TRACKR.dtrack = energy deposition of the jth deposition even
646 if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
650 for ( Int_t j=0;j<TRACKR.mtrack;j++) {
651 sum +=TRACKR.dtrack[j];
657 Int_t TFluka::TrackPid() const
659 // Return the id of the particle transported
660 // TRACKR.jtrack = identity number of the particle
661 return PDGFromId(TRACKR.jtrack);
664 Double_t TFluka::TrackCharge() const
666 // Return charge of the track currently transported
667 // PAPROP.ichrge = electric charge of the particle
668 // TRACKR.jtrack = identity number of the particle
669 return PAPROP.ichrge[TRACKR.jtrack+6];
672 Double_t TFluka::TrackMass() const
674 // PAPROP.am = particle mass in GeV
675 // TRACKR.jtrack = identity number of the particle
676 return PAPROP.am[TRACKR.jtrack+6];
679 Double_t TFluka::Etot() const
681 // TRACKR.etrack = total energy of the particle
682 return TRACKR.etrack;
688 Bool_t TFluka::IsNewTrack() const
691 // True if the track is not at the boundary of the current volume
692 // Not true in some cases in bxdraw - to be solved
696 Bool_t TFluka::IsTrackInside() const
698 // True if the track is not at the boundary of the current volume
699 // In Fluka a step is always inside one kind of material
700 // If the step would go behind the region of one material,
701 // it will be shortened to reach only the boundary.
702 // Therefore IsTrackInside() is always true.
703 // Not true in some cases in bxdraw - to be solved
707 Bool_t TFluka::IsTrackEntering() const
709 // True if this is the first step of the track in the current volume
710 // Boundary- (X) crossing
711 // Icode = 19: boundary crossing - call from Kaskad
712 // Icode = 29: boundary crossing - call from Emfsco
713 // Icode = 39: boundary crossing - call from Kasneu
714 // Icode = 49: boundary crossing - call from Kashea
715 // Icode = 59: boundary crossing - call from Kasoph
720 fIcode == 59) return 1;
724 Bool_t TFluka::IsTrackExiting() const
726 // True if this is the last step of the track in the current volume
727 // Boundary- (X) crossing
728 // Icode = 19: boundary crossing - call from Kaskad
729 // Icode = 29: boundary crossing - call from Emfsco
730 // Icode = 39: boundary crossing - call from Kasneu
731 // Icode = 49: boundary crossing - call from Kashea
732 // Icode = 59: boundary crossing - call from Kasoph
737 fIcode == 59) return 1;
741 Bool_t TFluka::IsTrackOut() const
743 // True if the track is out of the setup
745 // Icode = 14: escape - call from Kaskad
746 // Icode = 23: escape - call from Emfsco
747 // Icode = 32: escape - call from Kasneu
748 // Icode = 40: escape - call from Kashea
749 // Icode = 51: escape - call from Kasoph
754 fIcode == 51) return 1;
758 Bool_t TFluka::IsTrackDisappeared() const
760 // means all inelastic interactions and decays
761 // fIcode from usdraw
762 if (fIcode == 101 || // inelastic interaction
763 fIcode == 102 || // particle decay
764 fIcode == 214 || // in-flight annihilation
765 fIcode == 215 || // annihilation at rest
766 fIcode == 217 || // pair production
767 fIcode == 221) return 1;
771 Bool_t TFluka::IsTrackStop() const
773 // True if the track energy has fallen below the threshold
774 // means stopped by signal or below energy threshold
775 // Icode = 12: stopping particle - call from Kaskad
776 // Icode = 15: time kill - call from Kaskad
777 // Icode = 21: below threshold, iarg=1 - call from Emfsco
778 // Icode = 22: below threshold, iarg=2 - call from Emfsco
779 // Icode = 24: time kill - call from Emfsco
780 // Icode = 31: below threshold - call from Kasneu
781 // Icode = 33: time kill - call from Kasneu
782 // Icode = 41: time kill - call from Kashea
783 // Icode = 52: time kill - call from Kasoph
792 fIcode == 52) return 1;
796 Bool_t TFluka::IsTrackAlive() const
798 // means not disappeared or not out
799 if (IsTrackDisappeared() || IsTrackOut() ) return 0;
807 Int_t TFluka::NSecondaries() const
808 // Number of secondary particles generated in the current step
809 // FINUC.np = number of secondaries except light and heavy ions
810 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
812 return FINUC.np + FHEAVY.npheav;
815 void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
816 TLorentzVector& position, TLorentzVector& momentum)
818 if (isec >= 0 && isec < FINUC.np) {
819 // more fine condition depending on icode
837 particleId = PDGFromId(FINUC.kpart[isec]);
838 position.SetX(fXsco);
839 position.SetY(fYsco);
840 position.SetZ(fZsco);
841 position.SetT(TRACKR.atrack);
842 // position.SetT(TRACKR.atrack+FINUC.agesec[isec]); //not yet implem.
843 momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
844 momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
845 momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
846 momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
848 if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
849 Int_t jsec = isec - FINUC.np;
850 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
851 position.SetX(fXsco);
852 position.SetY(fYsco);
853 position.SetZ(fZsco);
854 position.SetT(TRACKR.atrack);
855 // position.SetT(TRACKR.atrack+FHEAVY.agheav[jsec]); //not yet implem.
856 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
857 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
858 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
859 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
860 momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
861 else if (FHEAVY.tkheav[jsec] > 6)
862 momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
866 TMCProcess TFluka::ProdProcess(Int_t isec) const
867 // Name of the process that has produced the secondary particles
868 // in the current step
870 const TMCProcess kIpNoProc = kPNoProcess;
871 const TMCProcess kIpPDecay = kPDecay;
872 const TMCProcess kIpPPair = kPPair;
873 //const TMCProcess kIpPPairFromPhoton = kPPairFromPhoton;
874 //const TMCProcess kIpPPairFromVirtualPhoton = kPPairFromVirtualPhoton;
875 const TMCProcess kIpPCompton = kPCompton;
876 const TMCProcess kIpPPhotoelectric = kPPhotoelectric;
877 const TMCProcess kIpPBrem = kPBrem;
878 //const TMCProcess kIpPBremFromHeavy = kPBremFromHeavy;
879 //const TMCProcess kIpPBremFromElectronOrPositron = kPBremFromElectronOrPositron;
880 const TMCProcess kIpPDeltaRay = kPDeltaRay;
881 //const TMCProcess kIpPMoller = kPMoller;
882 //const TMCProcess kIpPBhabha = kPBhabha;
883 const TMCProcess kIpPAnnihilation = kPAnnihilation;
884 //const TMCProcess kIpPAnnihilInFlight = kPAnnihilInFlight;
885 //const TMCProcess kIpPAnnihilAtRest = kPAnnihilAtRest;
886 const TMCProcess kIpPHadronic = kPHadronic;
887 const TMCProcess kIpPMuonNuclear = kPMuonNuclear;
888 const TMCProcess kIpPPhotoFission = kPPhotoFission;
889 const TMCProcess kIpPRayleigh = kPRayleigh;
890 // const TMCProcess kIpPCerenkov = kPCerenkov;
891 // const TMCProcess kIpPSynchrotron = kPSynchrotron;
893 Int_t mugamma = TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11;
894 if (fIcode == 102) return kIpPDecay;
895 else if (fIcode == 104 || fIcode == 217) return kIpPPair;
896 //else if (fIcode == 104) return kIpPairFromPhoton;
897 //else if (fIcode == 217) return kIpPPairFromVirtualPhoton;
898 else if (fIcode == 219) return kIpPCompton;
899 else if (fIcode == 221) return kIpPPhotoelectric;
900 else if (fIcode == 105 || fIcode == 208) return kIpPBrem;
901 //else if (fIcode == 105) return kIpPBremFromHeavy;
902 //else if (fIcode == 208) return kPBremFromElectronOrPositron;
903 else if (fIcode == 103 || fIcode == 400) return kIpPDeltaRay;
904 else if (fIcode == 210 || fIcode == 212) return kIpPDeltaRay;
905 //else if (fIcode == 210) return kIpPMoller;
906 //else if (fIcode == 212) return kIpPBhabha;
907 else if (fIcode == 214 || fIcode == 215) return kIpPAnnihilation;
908 //else if (fIcode == 214) return kIpPAnnihilInFlight;
909 //else if (fIcode == 215) return kIpPAnnihilAtRest;
910 else if (fIcode == 101) return kIpPHadronic;
911 else if (fIcode == 101) {
912 if (!mugamma) return kIpPHadronic;
913 else if (TRACKR.jtrack == 7) return kIpPPhotoFission;
914 else return kIpPMuonNuclear;
916 else if (fIcode == 225) return kIpPRayleigh;
917 // Fluka codes 100, 300 and 400 still to be investigasted
918 else return kIpNoProc;
921 //Int_t StepProcesses(TArrayI &proc) const
922 // Return processes active in the current step
924 //ck = total energy of the particl ????????????????
928 Int_t TFluka::VolId2Mate(Int_t id) const
931 // Returns the material number for a given volume ID
933 printf("VolId2Mate %d %d\n", id, fMediaByRegion[id]);
934 return fMediaByRegion[id-1];
937 const char* TFluka::VolName(Int_t id) const
940 // Returns the volume name for a given volume ID
942 FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[id-1]);
943 const char* name = vol->GetName();
944 printf("VolName %d %s \n", id, name);
948 Int_t TFluka::VolId(const Text_t* volName) const
951 // Converts from volume name to volume ID.
952 // Time consuming. (Only used during set-up)
953 // Could be replaced by hash-table
957 for (i = 0; i < fNVolumes; i++)
959 FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]);
960 TString name = vol->GetName();
961 strcpy(tmp, name.Data());
963 if (!strcmp(tmp, volName)) break;
971 Int_t TFluka::CurrentVolID(Int_t& copyNo) const
974 // Return the logical id and copy number corresponding to the current fluka region
976 int ir = fCurrentFlukaRegion;
977 int id = (FGeometryInit::GetInstance())->CurrentVolID(ir, copyNo);
978 printf("CurrentVolID: %d %d %d \n", ir, id, copyNo);
983 Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
986 // Return the logical id and copy number of off'th mother
987 // corresponding to the current fluka region
990 return CurrentVolID(copyNo);
992 int ir = fCurrentFlukaRegion;
993 int id = (FGeometryInit::GetInstance())->CurrentVolOffID(ir, off, copyNo);
995 printf("CurrentVolOffID: %d %d %d \n", ir, id, copyNo);
997 printf("CurrentVolOffID: Warning Mother not found !!!\n");
1002 const char* TFluka::CurrentVolName() const
1005 // Return the current volume name
1008 Int_t id = TFluka::CurrentVolID(copy);
1009 const char* name = TFluka::VolName(id);
1010 printf("CurrentVolumeName: %d %s \n", fCurrentFlukaRegion, name);
1014 const char* TFluka::CurrentVolOffName(Int_t off) const
1017 // Return the volume name of the off'th mother of the current volume
1020 Int_t id = TFluka::CurrentVolOffID(off, copy);
1021 const char* name = TFluka::VolName(id);
1022 printf("CurrentVolumeOffName: %d %s \n", fCurrentFlukaRegion, name);
1026 Int_t TFluka::CurrentMaterial(Float_t &a, Float_t &z,
1027 Float_t &dens, Float_t &radl, Float_t &absl) const
1030 // Return the current medium number
1033 Int_t id = TFluka::CurrentVolID(copy);
1034 Int_t med = TFluka::VolId2Mate(id);
1035 printf("CurrentMaterial: %d %d \n", fCurrentFlukaRegion, med);
1040 // ===============================================================
1041 void TFluka::FutoTest()
1043 Int_t icode, mreg, newreg, particleId;
1045 Double_t rull, xsco, ysco, zsco;
1046 TLorentzVector position, momentum;
1049 cout << " icode=" << icode << endl;
1051 cout << "TLorentzVector positionX=" << position.X()
1052 << "positionY=" << position.Y()
1053 << "positionZ=" << position.Z()
1054 << "timeT=" << position.T() << endl;
1055 cout << "TLorentzVector momentumX=" << momentum.X()
1056 << "momentumY=" << momentum.Y()
1057 << "momentumZ=" << momentum.Z()
1058 << "energyE=" << momentum.E() << endl;
1059 cout << "TrackPid=" << TrackPid() << endl;
1063 else if (icode > 0 && icode <= 5) {
1066 // medium = GetMedium();
1067 cout << " icode=" << icode
1069 // << " medium=" << medium
1071 TrackPosition(position);
1072 TrackMomentum(momentum);
1073 cout << "TLorentzVector positionX=" << position.X()
1074 << "positionY=" << position.Y()
1075 << "positionZ=" << position.Z()
1076 << "timeT=" << position.T() << endl;
1077 cout << "TLorentzVector momentumX=" << momentum.X()
1078 << "momentumY=" << momentum.Y()
1079 << "momentumZ=" << momentum.Z()
1080 << "energyE=" << momentum.E() << endl;
1081 cout << "TrackStep=" << TrackStep() << endl;
1082 cout << "TrackLength=" << TrackLength() << endl;
1083 cout << "TrackTime=" << TrackTime() << endl;
1084 cout << "Edep=" << Edep() << endl;
1085 cout << "TrackPid=" << TrackPid() << endl;
1086 cout << "TrackCharge=" << TrackCharge() << endl;
1087 cout << "TrackMass=" << TrackMass() << endl;
1088 cout << "Etot=" << Etot() << endl;
1089 cout << "IsNewTrack=" << IsNewTrack() << endl;
1090 cout << "IsTrackInside=" << IsTrackInside() << endl;
1091 cout << "IsTrackEntering=" << IsTrackEntering() << endl;
1092 cout << "IsTrackExiting=" << IsTrackExiting() << endl;
1093 cout << "IsTrackOut=" << IsTrackOut() << endl;
1094 cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
1095 cout << "IsTrackAlive=" << IsTrackAlive() << endl;
1098 else if((icode >= 10 && icode <= 15) ||
1099 (icode >= 20 && icode <= 24) ||
1100 (icode >= 30 && icode <= 33) ||
1101 (icode >= 40 && icode <= 41) ||
1102 (icode >= 50 && icode <= 52)) {
1105 // medium = GetMedium();
1110 cout << " icode=" << icode
1112 // << " medium=" << medium
1116 << " zsco=" << zsco << endl;
1117 TrackPosition(position);
1118 TrackMomentum(momentum);
1119 cout << "Edep=" << Edep() << endl;
1120 cout << "Etot=" << Etot() << endl;
1121 cout << "TrackPid=" << TrackPid() << endl;
1122 cout << "TrackCharge=" << TrackCharge() << endl;
1123 cout << "TrackMass=" << TrackMass() << endl;
1124 cout << "IsTrackOut=" << IsTrackOut() << endl;
1125 cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
1126 cout << "IsTrackStop=" << IsTrackStop() << endl;
1127 cout << "IsTrackAlive=" << IsTrackAlive() << endl;
1130 else if((icode >= 100 && icode <= 105) ||
1134 (icode >= 214 && icode <= 215) ||
1143 // medium = GetMedium();
1147 cout << " icode=" << icode
1149 // << " medium=" << medium
1152 << " zsco=" << zsco << endl;
1153 cout << "TrackPid=" << TrackPid() << endl;
1154 cout << "NSecondaries=" << NSecondaries() << endl;
1155 for (Int_t isec=0; isec< NSecondaries(); isec++) {
1156 TFluka::GetSecondary(isec, particleId, position, momentum);
1157 cout << "TLorentzVector positionX=" << position.X()
1158 << "positionY=" << position.Y()
1159 << "positionZ=" << position.Z()
1160 << "timeT=" << position.T() << endl;
1161 cout << "TLorentzVector momentumX=" << momentum.X()
1162 << "momentumY=" << momentum.Y()
1163 << "momentumZ=" << momentum.Z()
1164 << "energyE=" << momentum.E() << endl;
1165 cout << "TrackPid=" << particleId << endl;
1170 else if((icode == 19) ||
1176 // medium = GetMedium();
1177 newreg = GetNewreg();
1181 cout << " icode=" << icode
1183 // << " medium=" << medium
1184 << " newreg=" << newreg
1187 << " zsco=" << zsco << endl;
1190 // ====================================================================
1195 } // end of FutoTest