/*
$Log$
+Revision 1.12 2003/01/31 12:18:53 morsch
+Corrected indices. (E. Futo)
+
+Revision 1.9 2002/12/06 12:41:29 morsch
+Mess from last merge cleaned up.
+
+Revision 1.8 2002/12/06 12:28:44 morsch
+Region to media mapping corrected and improved.
+
+Revision 1.7 2002/12/06 12:21:32 morsch
+User stepping methods added (E. Futo)
+
Revision 1.6 2002/11/21 18:40:06 iglez2
Media handling added
#include <Riostream.h>
+#include "TClonesArray.h"
#include "TFluka.h"
#include "TCallf77.h" //For the fortran calls
#include "Fdblprc.h" //(DBLPRC) fluka common
#include "Fpaprop.h" //(PAPROP) fluka common
#include "Fpart.h" //(PART) fluka common
#include "Ftrackr.h" //(TRACKR) fluka common
+#include "Fpaprop.h" //(PAPROP) fluka common
#include "Ffheavy.h" //(FHEAVY) fluka common
#include "TVirtualMC.h"
#include "FGeometryInit.hh"
#include "TLorentzVector.h"
+#include "FlukaVolume.h"
// Fluka methods that may be needed.
#ifndef WIN32
if (fVerbosityLevel >=3)
cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
+
+ fVolumeMediaMap = new TClonesArray("FlukaVolume",1000);
+ fNVolumes = 0;
+ fMediaByRegion = 0;
}
TFluka::~TFluka() {
cout << "==> TFluka::~TFluka() destructor called." << endl;
delete fGeometryManager;
+ fVolumeMediaMap->Delete();
+ delete fVolumeMediaMap;
+
if (fVerbosityLevel >=3)
cout << "<== TFluka::~TFluka() destructor called." << endl;
if (fVerbosityLevel >=3)
cout << "<== TFluka::Init() called." << endl;
+ FinishGeometry();
+
}
void TFluka::FinishGeometry() {
+//
+// Build-up table with region to medium correspondance
+//
+ char tmp[5];
+
if (fVerbosityLevel >=3)
cout << "==> TFluka::FinishGeometry() called." << endl;
- fGeometryManager->Ggclos();
+// fGeometryManager->Ggclos();
- FGeometryInit* flugg = FGeometryInit::GetInstance();
- map<TString, Int_t, less<TString> >::iterator i;
- for (fVolumeMediaMap.begin(); i != fVolumeMediaMap.end(); i++) {
- TString volName = (*i).first;
- Int_t media = (*i).second;
- Int_t region = flugg->GetRegionFromName(volName);
- fMediaByRegion[region] = media;
+ FGeometryInit* flugg = FGeometryInit::GetInstance();
+
+ fMediaByRegion = new Int_t[fNVolumes+2];
+ for (Int_t i = 0; i < fNVolumes; i++)
+ {
+ FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]);
+ TString volName = vol->GetName();
+ Int_t media = vol->GetMedium();
+ printf("Finish Geometry: volName, media %d %s %d \n", i, volName.Data(), media);
+ strcpy(tmp, volName.Data());
+ tmp[4] = '\0';
+ flugg->SetMediumFromName(tmp, media, i+1);
+ fMediaByRegion[i] = media;
}
+
+ flugg->BuildMediaMap();
if (fVerbosityLevel >=3)
cout << "<== TFluka::FinishGeometry() called." << endl;
void TFluka::ProcessEvent() {
if (fVerbosityLevel >=3)
cout << "==> TFluka::ProcessEvent() called." << endl;
-
+ fApplication->GeneratePrimaries();
+ EPISOR.lsouit = true;
+ flukam(1);
if (fVerbosityLevel >=3)
cout << "<== TFluka::ProcessEvent() called." << endl;
}
cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
cout << "\t* Calling flukam again..." << endl;
}
- fApplication->GeneratePrimaries();
- EPISOR.lsouit = true;
- flukam(1);
-
+ fApplication->InitGeometry();
+ fApplication->BeginEvent();
+ ProcessEvent();
+ fApplication->FinishEvent();
if (fVerbosityLevel >=3)
cout << "<== TFluka::ProcessRun(" << nevent << ") called."
<< endl;
+
}
//_____________________________________________________________________________
Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
Float_t *upar, Int_t np) {
//
- fVolumeMediaMap[TString(name)] = nmed;
- return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
+// fVolumeMediaMap[TString(name)] = nmed;
+ printf("TFluka::Gsvolu() name = %s, nmed = %d\n", name, nmed);
+
+ TClonesArray &lvols = *fVolumeMediaMap;
+ new(lvols[fNVolumes++])
+ FlukaVolume(name, nmed);
+ return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
}
Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
Double_t *upar, Int_t np) {
//
- return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
+ TClonesArray &lvols = *fVolumeMediaMap;
+ new(lvols[fNVolumes++])
+ FlukaVolume(name, nmed);
+
+ return fGeometryManager->Gsvolu(name, shape, nmed, upar, np);
}
void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
Int_t iaxis) {
//
- fGeometryManager->Gsdvn(name, mother, ndiv, iaxis);
+// The medium of the daughter is the one of the mother
+ Int_t volid = TFluka::VolId(mother);
+ Int_t med = TFluka::VolId2Mate(volid);
+ TClonesArray &lvols = *fVolumeMediaMap;
+ new(lvols[fNVolumes++])
+ FlukaVolume(name, med);
+ fGeometryManager->Gsdvn(name, mother, ndiv, iaxis);
}
void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
Int_t iaxis, Double_t c0i, Int_t numed) {
//
- fGeometryManager->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
+ TClonesArray &lvols = *fVolumeMediaMap;
+ new(lvols[fNVolumes++])
+ FlukaVolume(name, numed);
+ fGeometryManager->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
}
void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
Int_t iaxis, Int_t numed, Int_t ndvmx) {
-//
- fGeometryManager->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
+//
+ TClonesArray &lvols = *fVolumeMediaMap;
+ new(lvols[fNVolumes++])
+ FlukaVolume(name, numed);
+ fGeometryManager->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
}
void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
//
- fGeometryManager->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
+ TClonesArray &lvols = *fVolumeMediaMap;
+ new(lvols[fNVolumes++])
+ FlukaVolume(name, numed);
+ fGeometryManager->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
}
void TFluka::Gsord(const char *name, Int_t iax) {
//_____________________________________________________________________________
// methods needed by the stepping
//____________________________________________________________________________
+
Int_t TFluka::GetMedium() const {
- return fMediaByRegion[fCurrentFlukaRegion];
+//
+// Get the medium number for the current fluka region
+//
+ FGeometryInit* flugg = FGeometryInit::GetInstance();
+ return flugg->GetMedium(fCurrentFlukaRegion);
}
// Return PDG code and pseudo ENDF code from Fluka code
//IPTOKP array goes from official to internal
+ if (id == 0) {
+ printf("PDGFromId: Error id = 0");
+ return -1;
+ }
+
Int_t intfluka = GetFlukaIPTOKP(id);
- //MPKDHA() goes from internal to PDG
+ if (intfluka == 0) {
+ printf("PDGFromId: Error intfluka = 0");
+ return -1;
+ }
+
+ //MPKDHA() goes from internal to PDG
return mpdgha(intfluka);
-
}
-
-
//_____________________________________________________________________________
// methods for step management
//____________________________________________________________________________
+//
+// set methods
+//
+void TFluka::SetMaxStep(Double_t)
+{
+// SetMaxStep is dummy procedure in TFluka !
+ cout << "SetMaxStep is dummy procedure in TFluka !" << endl;
+}
+
+void TFluka::SetMaxNStep(Int_t)
+{
+// SetMaxNStep is dummy procedure in TFluka !
+ cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
+}
+
+void TFluka::SetUserDecay(Int_t)
+{
+// SetUserDecay is dummy procedure in TFluka !
+ cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
+}
+
//
// dynamic properties
//
return;
}
else {
- Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack]*PAPROP.am[TRACKR.jtrack]);
+ Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
momentum.SetPx(p*TRACKR.cxtrck);
momentum.SetPy(p*TRACKR.cytrck);
momentum.SetPz(p*TRACKR.cztrck);
Double_t TFluka::TrackLength() const
{
-// It is wrong
-// should be the sum of all steps starting from the beginning of the track
+// Still wrong !!!
+// This is the sum of substeps !!!
+// TRACKR.ctrack = total curved path of the current step
+// Sum of the substeps is identical to TRACKR.ctrack if the is no mag. field
+// The sum of all step length starting from the beginning of the track
// for the time being returns only the length in centimeters of the current step
- return TRACKR.ctrack;
+ Double_t sum = 0;
+ for ( Int_t j=0;j<TRACKR.ntrack;j++) {
+ sum +=TRACKR.ttrack[j];
+ }
+ return sum;
}
Double_t TFluka::TrackTime() const
{
// Return charge of the track currently transported
// PAPROP.ichrge = electric charge of the particle
+// TRACKR.jtrack = identity number of the particle
return PAPROP.ichrge[TRACKR.jtrack+6];
}
Double_t TFluka::TrackMass() const
{
// PAPROP.am = particle mass in GeV
+// TRACKR.jtrack = identity number of the particle
return PAPROP.am[TRACKR.jtrack+6];
}
Int_t TFluka::NSecondaries() const
// Number of secondary particles generated in the current step
-// fIcode = 100 = elastic interaction
-// fIcode = 101 = inelastic interaction
-// fIcode = 102 = particle decay
-// fIcode = 103 = delta ray
-// fIcode = 104 = pair production
-// fIcode = 105 = bremsstrahlung
+// FINUC.np = number of secondaries except light and heavy ions
+// FHEAVY.npheav = number of secondaries for light and heavy secondary ions
{
- if (fIcode >= 100 && fIcode <= 105)
- return FINUC.np + FHEAVY.npheav;
- else
- return -1;
+ return FINUC.np + FHEAVY.npheav;
}
void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
TLorentzVector& position, TLorentzVector& momentum)
-//
-// fIcode = 100 = elastic interaction
-// fIcode = 101 = inelastic interaction
-// fIcode = 102 = particle decay
-// fIcode = 103 = delta ray
-// fIcode = 104 = pair production
-// fIcode = 105 = bremsstrahlung
-{
- if (fIcode >= 100 && fIcode <= 105) {
- if (isec >= 0 && isec < FINUC.np) {
- particleId = PDGFromId(FINUC.kpart[isec]);
- position.SetX(fXsco);
- position.SetY(fYsco);
- position.SetZ(fZsco);
- position.SetT(FINUC.agesec[isec]);
- momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
- momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
- momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
- momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
- }
- if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
- Int_t jsec = isec - FINUC.np;
- particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
- position.SetX(fXsco);
- position.SetY(fYsco);
- position.SetZ(fZsco);
- position.SetT(FHEAVY.agheav[jsec]);
- momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
- momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
- momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
- if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
- momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
- else if (FHEAVY.tkheav[jsec] > 6)
- momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
- }
+{
+ if (isec >= 0 && isec < FINUC.np) {
+ // more fine condition depending on icode
+ // icode = 100 ?
+ // icode = 101 OK
+ // icode = 102 OK
+ // icode = 103 ?
+ // icode = 104 ?
+ // icode = 105 ?
+ // icode = 208 ?
+ // icode = 210 ?
+ // icode = 212 ?
+ // icode = 214 OK
+ // icode = 215 OK
+ // icode = 219 ?
+ // icode = 221 OK
+ // icode = 225 ?
+ // icode = 300 ?
+ // icode = 400 ?
+
+ particleId = PDGFromId(FINUC.kpart[isec]);
+ position.SetX(fXsco);
+ position.SetY(fYsco);
+ position.SetZ(fZsco);
+ position.SetT(TRACKR.atrack);
+// position.SetT(TRACKR.atrack+FINUC.agesec[isec]); //not yet implem.
+ momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
+ momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
+ momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
+ momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
}
+ if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
+ Int_t jsec = isec - FINUC.np;
+ particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
+ position.SetX(fXsco);
+ position.SetY(fYsco);
+ position.SetZ(fZsco);
+ position.SetT(TRACKR.atrack);
+// position.SetT(TRACKR.atrack+FHEAVY.agheav[jsec]); //not yet implem.
+ momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
+ momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
+ momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
+ if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
+ momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
+ else if (FHEAVY.tkheav[jsec] > 6)
+ momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
+ }
}
-//TMCProcess ProdProcess(Int_t isec) const
+TMCProcess TFluka::ProdProcess(Int_t isec) const
// Name of the process that has produced the secondary particles
// in the current step
-//{
-// will come from FINUC when called from USDRAW
-//}
+{
+ const TMCProcess kIpNoProc = kPNoProcess;
+ const TMCProcess kIpPDecay = kPDecay;
+ const TMCProcess kIpPPair = kPPair;
+//const TMCProcess kIpPPairFromPhoton = kPPairFromPhoton;
+//const TMCProcess kIpPPairFromVirtualPhoton = kPPairFromVirtualPhoton;
+ const TMCProcess kIpPCompton = kPCompton;
+ const TMCProcess kIpPPhotoelectric = kPPhotoelectric;
+ const TMCProcess kIpPBrem = kPBrem;
+//const TMCProcess kIpPBremFromHeavy = kPBremFromHeavy;
+//const TMCProcess kIpPBremFromElectronOrPositron = kPBremFromElectronOrPositron;
+ const TMCProcess kIpPDeltaRay = kPDeltaRay;
+//const TMCProcess kIpPMoller = kPMoller;
+//const TMCProcess kIpPBhabha = kPBhabha;
+ const TMCProcess kIpPAnnihilation = kPAnnihilation;
+//const TMCProcess kIpPAnnihilInFlight = kPAnnihilInFlight;
+//const TMCProcess kIpPAnnihilAtRest = kPAnnihilAtRest;
+ const TMCProcess kIpPHadronic = kPHadronic;
+ const TMCProcess kIpPMuonNuclear = kPMuonNuclear;
+ const TMCProcess kIpPPhotoFission = kPPhotoFission;
+ const TMCProcess kIpPRayleigh = kPRayleigh;
+// const TMCProcess kIpPCerenkov = kPCerenkov;
+// const TMCProcess kIpPSynchrotron = kPSynchrotron;
+
+ Int_t mugamma = TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11;
+ if (fIcode == 102) return kIpPDecay;
+ else if (fIcode == 104 || fIcode == 217) return kIpPPair;
+//else if (fIcode == 104) return kIpPairFromPhoton;
+//else if (fIcode == 217) return kIpPPairFromVirtualPhoton;
+ else if (fIcode == 219) return kIpPCompton;
+ else if (fIcode == 221) return kIpPPhotoelectric;
+ else if (fIcode == 105 || fIcode == 208) return kIpPBrem;
+//else if (fIcode == 105) return kIpPBremFromHeavy;
+//else if (fIcode == 208) return kPBremFromElectronOrPositron;
+ else if (fIcode == 103 || fIcode == 400) return kIpPDeltaRay;
+ else if (fIcode == 210 || fIcode == 212) return kIpPDeltaRay;
+//else if (fIcode == 210) return kIpPMoller;
+//else if (fIcode == 212) return kIpPBhabha;
+ else if (fIcode == 214 || fIcode == 215) return kIpPAnnihilation;
+//else if (fIcode == 214) return kIpPAnnihilInFlight;
+//else if (fIcode == 215) return kIpPAnnihilAtRest;
+ else if (fIcode == 101) return kIpPHadronic;
+ else if (fIcode == 101) {
+ if (!mugamma) return kIpPHadronic;
+ else if (TRACKR.jtrack == 7) return kIpPPhotoFission;
+ else return kIpPMuonNuclear;
+ }
+ else if (fIcode == 225) return kIpPRayleigh;
+// Fluka codes 100, 300 and 400 still to be investigasted
+ else return kIpNoProc;
+}
//Int_t StepProcesses(TArrayI &proc) const
// Return processes active in the current step
//}
+Int_t TFluka::VolId2Mate(Int_t id) const
+{
+//
+// Returns the material number for a given volume ID
+//
+ printf("VolId2Mate %d %d\n", id, fMediaByRegion[id]);
+ return fMediaByRegion[id-1];
+}
+
+const char* TFluka::VolName(Int_t id) const
+{
+//
+// Returns the volume name for a given volume ID
+//
+ FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[id-1]);
+ const char* name = vol->GetName();
+ printf("VolName %d %s \n", id, name);
+ return name;
+}
+
+Int_t TFluka::VolId(const Text_t* volName) const
+{
+//
+// Converts from volume name to volume ID.
+// Time consuming. (Only used during set-up)
+// Could be replaced by hash-table
+//
+ char tmp[5];
+ Int_t i =0;
+ for (i = 0; i < fNVolumes; i++)
+ {
+ FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]);
+ TString name = vol->GetName();
+ strcpy(tmp, name.Data());
+ tmp[4] = '\0';
+ if (!strcmp(tmp, volName)) break;
+ }
+ i++;
+
+ return i;
+}
+
+
+Int_t TFluka::CurrentVolID(Int_t& copyNo) const
+{
+//
+// Return the logical id and copy number corresponding to the current fluka region
+//
+ int ir = fCurrentFlukaRegion;
+ int id = (FGeometryInit::GetInstance())->CurrentVolID(ir, copyNo);
+ printf("CurrentVolID: %d %d %d \n", ir, id, copyNo);
+ return id;
+
+}
+
+Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
+{
+//
+// Return the logical id and copy number of off'th mother
+// corresponding to the current fluka region
+//
+ if (off == 0)
+ return CurrentVolID(copyNo);
+
+ int ir = fCurrentFlukaRegion;
+ int id = (FGeometryInit::GetInstance())->CurrentVolOffID(ir, off, copyNo);
+
+ printf("CurrentVolOffID: %d %d %d \n", ir, id, copyNo);
+ if (id == -1)
+ printf("CurrentVolOffID: Warning Mother not found !!!\n");
+ return id;
+}
+
+
+const char* TFluka::CurrentVolName() const
+{
+//
+// Return the current volume name
+//
+ Int_t copy;
+ Int_t id = TFluka::CurrentVolID(copy);
+ const char* name = TFluka::VolName(id);
+ printf("CurrentVolumeName: %d %s \n", fCurrentFlukaRegion, name);
+ return name;
+}
+
+const char* TFluka::CurrentVolOffName(Int_t off) const
+{
+//
+// Return the volume name of the off'th mother of the current volume
+//
+ Int_t copy;
+ Int_t id = TFluka::CurrentVolOffID(off, copy);
+ const char* name = TFluka::VolName(id);
+ printf("CurrentVolumeOffName: %d %s \n", fCurrentFlukaRegion, name);
+ return name;
+}
+
+Int_t TFluka::CurrentMaterial(Float_t &a, Float_t &z,
+ Float_t &dens, Float_t &radl, Float_t &absl) const
+{
+//
+// Return the current medium number
+//
+ Int_t copy;
+ Int_t id = TFluka::CurrentVolID(copy);
+ Int_t med = TFluka::VolId2Mate(id);
+ printf("CurrentMaterial: %d %d \n", fCurrentFlukaRegion, med);
+ return med;
+}
+
+
// ===============================================================
void TFluka::FutoTest()
{
cout << "TrackPid=" << TrackPid() << endl;
cout << "NSecondaries=" << NSecondaries() << endl;
for (Int_t isec=0; isec< NSecondaries(); isec++) {
-//void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
-// TLorentzVector& position, TLorentzVector& momentum)
TFluka::GetSecondary(isec, particleId, position, momentum);
cout << "TLorentzVector positionX=" << position.X()
<< "positionY=" << position.Y()
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff