//
#include <Riostream.h>
+#include <TList.h>
#include "TFluka.h"
#include "TFlukaCodes.h"
fTrackIsExiting(kFALSE),
fTrackIsNew(kFALSE),
fFieldFlag(kTRUE),
- fGeneratePemf(kFALSE),
fDummyBoundary(kFALSE),
fStopped(kFALSE),
fStopEvent(kFALSE),
fStopRun(kFALSE),
+ fPrimaryElectronIndex(-1),
fMaterials(0),
fNVolumes(0),
fCurrentFlukaRegion(-1),
//
// Default constructor
//
+ for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
}
//______________________________________________________________________________
fTrackIsExiting(kFALSE),
fTrackIsNew(kFALSE),
fFieldFlag(kTRUE),
- fGeneratePemf(kFALSE),
fDummyBoundary(kFALSE),
fStopped(kFALSE),
fStopEvent(kFALSE),
fStopRun(kFALSE),
+ fPrimaryElectronIndex(-1),
fMaterials(0),
fNVolumes(0),
fCurrentFlukaRegion(-1),
fUserScore(new TObjArray(100))
{
// create geometry interface
+ for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
+
if (fVerbosityLevel >=3)
cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
SetCoreInputFileName();
cout << "<== TFluka::~TFluka() destructor called." << endl;
if (fMaterials) delete [] fMaterials;
- delete fGeom;
- delete fMCGeo;
+// delete fGeom;
+// delete fMCGeo;
if (fUserConfig) {
fUserConfig->Delete();
}
fApplication->InitGeometry();
-
+ fApplication->ConstructOpGeometry();
//
// Add ions to PDG Data base
//
AddParticlesToPdgDataBase();
+ //
}
//
if (fVerbosityLevel >=3) {
cout << "==> TFluka::FinishGeometry() called." << endl;
- printf("----FinishGeometry - nothing to do with TGeo\n");
+ printf("----FinishGeometry - applying misalignment if any\n");
cout << "<== TFluka::FinishGeometry() called." << endl;
}
+ TVirtualMCApplication::Instance()->MisalignGeometry();
}
//______________________________________________________________________________
}
}
- //
- // At this stage we have the information on materials and cuts available.
- // Now create the pemf file
-
- if (fGeneratePemf) fGeom->CreatePemfFile();
-
- //
+
// Prepare input file with the current physics settings
InitPhysics();
//
// Run steering
//
-
+
if (fVerbosityLevel >=3)
cout << "==> TFluka::ProcessRun(" << nevent << ") called."
<< endl;
// In this case, WMAT in output is changed to relative
// weigths.
//
+ printf("Mixture %5d %10s %5d \n", kmat, name, nlmat);
+
Int_t i,j;
if (nlmat < 0) {
nlmat = - nlmat;
delete [] wmatnew;
return;
}
- gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
+ printf("Mixture (2) %5d %10s %5d \n", kmat, name, nlmat);
+ gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
}
//______________________________________________________________________________
//
//
// Create object holding Cerenkov properties
-//
+//
+
TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
//
// Pass object to medium
//______________________________________________________________________________
-void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/,
- Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
+void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
+ Double_t *absco, Double_t *effic, Double_t *rindex) {
+//
+// Set Cerenkov properties for medium itmed
//
-// Double_t version not implemented
+// npckov: number of sampling points
+// ppckov: energy values
+// absco: absorption length
+// effic: quantum efficiency
+// rindex: refraction index
+//
+
+//
+// Double_t version
+ Float_t* fppckov = CreateFloatArray(ppckov, npckov);
+ Float_t* fabsco = CreateFloatArray(absco, npckov);
+ Float_t* feffic = CreateFloatArray(effic, npckov);
+ Float_t* frindex = CreateFloatArray(rindex, npckov);
+
+ SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex);
+
+ delete [] fppckov;
+ delete [] fabsco;
+ delete [] feffic;
+ delete [] frindex;
}
-void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t* /*ppckov*/,
- Double_t* /*absco*/, Double_t* /*effic*/, Double_t* /*rindex*/, Double_t* /*rfl*/) {
+void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t* ppckov,
+ Double_t* absco, Double_t* effic, Double_t* rindex, Double_t* rfl) {
+//
+// Set Cerenkov properties for medium itmed
//
-// // Double_t version not implemented
+// npckov: number of sampling points
+// ppckov: energy values
+// absco: absorption length
+// effic: quantum efficiency
+// rindex: refraction index
+// rfl: reflectivity for boundary to medium itmed
+//
+
+//
+// // Double_t version
+ Float_t* fppckov = CreateFloatArray(ppckov, npckov);
+ Float_t* fabsco = CreateFloatArray(absco, npckov);
+ Float_t* feffic = CreateFloatArray(effic, npckov);
+ Float_t* frindex = CreateFloatArray(rindex, npckov);
+ Float_t* frfl = CreateFloatArray(rfl, npckov);
+
+ SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex, frfl);
+
+ delete [] fppckov;
+ delete [] fabsco;
+ delete [] feffic;
+ delete [] frindex;
+ delete [] frfl;
}
// Euclid
//
// Get the medium number for the current fluka region
//
- return fGeom->GetMedium(); // this I need to check due to remapping !!!
+ if (gGeoManager->IsOutside()) {
+ return (-1);
+ } else {
+ return (fGeom->GetMedium()); // this I need to check due to remapping !!!
+ }
}
//____________________________________________________________________________
//
// Return PDG code and pseudo ENDF code from Fluka code
// Alpha He3 Triton Deuteron gen. ion opt. photon
- Int_t idSpecial[6] = {10020040, 10020030, 10010030, 10010020, 10000000, 50000050};
+ Int_t idSpecial[6] = {GetIonPdg(2,4), GetIonPdg(2, 3), GetIonPdg(1,3), GetIonPdg(1,2), GetIonPdg(0,0), 50000050};
// IPTOKP array goes from official to internal
if (id == kFLUKAoptical) {
// Error id
if (id == 0 || id < kFLUKAcodemin || id > kFLUKAcodemax) {
if (fVerbosityLevel >= 3)
- printf("PDGFromId: Error id = 0\n");
+ printf("PDGFromId: Error id = 0 %5d %5d\n", id, fCaller);
return -1;
}
// Good id
}
-void TFluka::SetUserScoring(const char* option, Int_t npr, char* outfile, Float_t* what)
+void TFluka::SetUserScoring(const char* option, const char* sdum, Int_t npr, char* outfile, Float_t* what)
{
//
// Adds a user scoring option to the list
//
- TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npr,outfile,what);
+ TFlukaScoringOption* opt = new TFlukaScoringOption(option, sdum, npr,outfile,what);
fUserScore->Add(opt);
}
//______________________________________________________________________________
-void TFluka::SetUserScoring(const char* option, Int_t npr, char* outfile, Float_t* what, const char* det1, const char* det2, const char* det3)
+void TFluka::SetUserScoring(const char* option, const char* sdum, Int_t npr, char* outfile, Float_t* what,
+ const char* det1, const char* det2, const char* det3)
{
//
// Adds a user scoring option to the list
//
- TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npr, outfile, what, det1, det2, det3);
+ TFlukaScoringOption* opt = new TFlukaScoringOption(option, sdum, npr, outfile, what, det1, det2, det3);
fUserScore->Add(opt);
}
// Process Fluka specific scoring options
//
TFlukaScoringOption::SetStaticInfo(pFlukaVmcInp, fGeom);
- Float_t loginp = 49.0;
+ Float_t loginp = -49.0;
Int_t inp = 0;
Int_t nscore = fUserScore->GetEntries();
// Initialisation needed for Cerenkov photon production and transport
TObjArray *matList = GetFlukaMaterials();
Int_t nmaterial = matList->GetEntriesFast();
- fMaterials = new Int_t[nmaterial+3];
+ fMaterials = new Int_t[nmaterial+25];
for (Int_t im = 0; im < nmaterial; im++)
{
if (caller == kENDRAW || caller == kUSDRAW ||
caller == kBXExiting || caller == kBXEntering ||
caller == kUSTCKV) {
- position.SetX(GetXsco());
- position.SetY(GetYsco());
- position.SetZ(GetZsco());
- position.SetT(TRACKR.atrack);
+ position.SetX(GetXsco());
+ position.SetY(GetYsco());
+ position.SetZ(GetZsco());
+ position.SetT(TRACKR.atrack);
}
- else if (caller == kMGDRAW) {
- position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
- position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
- position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
- position.SetT(TRACKR.atrack);
+ else if (caller == kMGDRAW) {
+ Int_t i = -1;
+ if ((i = fPrimaryElectronIndex) > -1) {
+ // Primary Electron Ionisation
+ Double_t x, y, z;
+ GetPrimaryElectronPosition(i, x, y, z);
+ position.SetX(x);
+ position.SetY(y);
+ position.SetZ(z);
+ position.SetT(TRACKR.atrack);
+ } else {
+ position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
+ position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
+ position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
+ position.SetT(TRACKR.atrack);
+ }
}
else if (caller == kSODRAW) {
- position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
- position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
- position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
- position.SetT(0);
+ Int_t ist = FLKSTK.npflka;
+ position.SetX(FLKSTK.xflk[ist]);
+ position.SetY(FLKSTK.yflk[ist]);
+ position.SetZ(FLKSTK.zflk[ist]);
+ position.SetT(FLKSTK.agestk[ist]);
} else if (caller == kMGResumedTrack) {
- position.SetX(TRACKR.spausr[0]);
- position.SetY(TRACKR.spausr[1]);
- position.SetZ(TRACKR.spausr[2]);
- position.SetT(TRACKR.spausr[3]);
+ position.SetX(TRACKR.spausr[0]);
+ position.SetY(TRACKR.spausr[1]);
+ position.SetZ(TRACKR.spausr[2]);
+ position.SetT(TRACKR.spausr[3]);
}
else
- Warning("TrackPosition","position not available");
+ Warning("TrackPosition","position not available");
}
//______________________________________________________________________________
if (caller == kENDRAW || caller == kUSDRAW ||
caller == kBXExiting || caller == kBXEntering ||
caller == kUSTCKV) {
- x = GetXsco();
- y = GetYsco();
- z = GetZsco();
+ x = GetXsco();
+ y = GetYsco();
+ z = GetZsco();
}
- else if (caller == kMGDRAW || caller == kSODRAW) {
- x = TRACKR.xtrack[TRACKR.ntrack];
- y = TRACKR.ytrack[TRACKR.ntrack];
- z = TRACKR.ztrack[TRACKR.ntrack];
+ else if (caller == kMGDRAW) {
+ Int_t i = -1;
+ if ((i = fPrimaryElectronIndex) > -1) {
+ GetPrimaryElectronPosition(i, x, y, z);
+ } else {
+ x = TRACKR.xtrack[TRACKR.ntrack];
+ y = TRACKR.ytrack[TRACKR.ntrack];
+ z = TRACKR.ztrack[TRACKR.ntrack];
+ }
+ }
+ else if (caller == kSODRAW) {
+ Int_t ist = FLKSTK.npflka;
+ x = FLKSTK.xflk[ist];
+ y = FLKSTK.yflk[ist];
+ z = FLKSTK.zflk[ist];
}
else if (caller == kMGResumedTrack) {
- x = TRACKR.spausr[0];
- y = TRACKR.spausr[1];
- z = TRACKR.spausr[2];
+ x = TRACKR.spausr[0];
+ y = TRACKR.spausr[1];
+ z = TRACKR.spausr[2];
}
else
- Warning("TrackPosition","position not available");
+ Warning("TrackPosition","position not available");
}
//______________________________________________________________________________
FlukaCallerCode_t caller = GetCaller();
FlukaProcessCode_t icode = GetIcode();
- if (caller != kEEDRAW && caller != kMGResumedTrack &&
+ if (caller != kEEDRAW &&
+ caller != kMGResumedTrack &&
+ caller != kSODRAW &&
+ caller != kUSDRAW &&
(caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
- if (TRACKR.ptrack >= 0) {
- momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
- momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
- momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
- momentum.SetE(TRACKR.etrack);
- return;
- }
- else {
- Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
- momentum.SetPx(p*TRACKR.cxtrck);
- momentum.SetPy(p*TRACKR.cytrck);
- momentum.SetPz(p*TRACKR.cztrck);
- momentum.SetE(TRACKR.etrack);
- return;
- }
+ if (TRACKR.ptrack >= 0) {
+ momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
+ momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
+ momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
+ momentum.SetE(TRACKR.etrack);
+ return;
+ }
+ else {
+ Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
+ momentum.SetPx(p*TRACKR.cxtrck);
+ momentum.SetPy(p*TRACKR.cytrck);
+ momentum.SetPz(p*TRACKR.cztrck);
+ momentum.SetE(TRACKR.etrack);
+ return;
+ }
} else if (caller == kMGResumedTrack) {
- momentum.SetPx(TRACKR.spausr[4]);
- momentum.SetPy(TRACKR.spausr[5]);
- momentum.SetPz(TRACKR.spausr[6]);
- momentum.SetE (TRACKR.spausr[7]);
- return;
+ momentum.SetPx(TRACKR.spausr[4]);
+ momentum.SetPy(TRACKR.spausr[5]);
+ momentum.SetPz(TRACKR.spausr[6]);
+ momentum.SetE (TRACKR.spausr[7]);
+ return;
} else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
momentum.SetPx(0.);
momentum.SetPy(0.);
momentum.SetPz(0.);
momentum.SetE(TrackMass());
+
+ } else if (caller == kSODRAW) {
+ Int_t ist = FLKSTK.npflka;
+ Double_t p = FLKSTK.pmoflk[ist];
+ Int_t ifl = FLKSTK.iloflk[ist];
+ Double_t m = PAPROP.am[ifl + 6];
+ Double_t e = TMath::Sqrt(p * p + m * m);
+ momentum.SetPx(p * FLKSTK.txflk[ist]);
+ momentum.SetPy(p * FLKSTK.tyflk[ist]);
+ momentum.SetPz(p * FLKSTK.tzflk[ist]);
+ momentum.SetE(e);
+ } else if (caller == kUSDRAW) {
+ if (icode == kEMFSCObrems ||
+ icode == kEMFSCOmoller ||
+ icode == kEMFSCObhabha ||
+ icode == kEMFSCOcompton )
+ {
+ momentum.SetPx(fPint[0]);
+ momentum.SetPy(fPint[1]);
+ momentum.SetPz(fPint[2]);
+ momentum.SetE(fPint[3]);
+ } else if (icode == kKASKADdray ||
+ icode == kKASKADbrems ||
+ icode == kKASKADpair) {
+ momentum.SetPx(GENSTK.plr[0] * GENSTK.cxr[0]);
+ momentum.SetPy(GENSTK.plr[0] * GENSTK.cyr[0]);
+ momentum.SetPz(GENSTK.plr[0] * GENSTK.czr[0]);
+ momentum.SetE (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);
+ } else {
+ Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack
+ - ParticleMassFPC(TRACKR.jtrack)
+ * ParticleMassFPC(TRACKR.jtrack));
+ momentum.SetPx(p*TRACKR.cxtrck);
+ momentum.SetPy(p*TRACKR.cytrck);
+ momentum.SetPz(p*TRACKR.cztrck);
+ momentum.SetE(TRACKR.etrack);
+ }
}
else
Warning("TrackMomentum","momentum not available");
// PAPROP.am[TRACKR.jtrack] = particle mass in gev
FlukaCallerCode_t caller = GetCaller();
FlukaProcessCode_t icode = GetIcode();
- if (caller != kEEDRAW && caller != kMGResumedTrack &&
+ if (caller != kEEDRAW &&
+ caller != kMGResumedTrack &&
+ caller != kSODRAW &&
+ caller != kUSDRAW &&
(caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
if (TRACKR.ptrack >= 0) {
px = TRACKR.ptrack*TRACKR.cxtrck;
py = 0.;
pz = 0.;
e = TrackMass();
+ } else if (caller == kSODRAW) {
+ Int_t ist = FLKSTK.npflka;
+ Double_t p = FLKSTK.pmoflk[ist];
+ Int_t ifl = FLKSTK.iloflk[ist];
+ Double_t m = PAPROP.am[ifl + 6];
+ e = TMath::Sqrt(p * p + m * m);
+ px = p * FLKSTK.txflk[ist];
+ py = p * FLKSTK.tyflk[ist];
+ pz = p * FLKSTK.tzflk[ist];
+ } else if (caller == kUSDRAW) {
+ if (icode == kEMFSCObrems ||
+ icode == kEMFSCOmoller ||
+ icode == kEMFSCObhabha ||
+ icode == kEMFSCOcompton )
+ {
+ px = fPint[0];
+ py = fPint[1];
+ pz = fPint[2];
+ e = fPint[3];
+ } else if (icode == kKASKADdray ||
+ icode == kKASKADbrems ||
+ icode == kKASKADpair) {
+ px = GENSTK.plr[0] * GENSTK.cxr[0];
+ py = GENSTK.plr[0] * GENSTK.cyr[0];
+ pz = GENSTK.plr[0] * GENSTK.czr[0];
+ e = GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6];
+ } else {
+ Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
+ px = p*TRACKR.cxtrck;
+ py = p*TRACKR.cytrck;
+ pz = p*TRACKR.cztrck;
+ e = TRACKR.etrack;
+ }
}
else
- Warning("TrackMomentum","momentum not available");
+ Warning("TrackMomentum","momentum not available");
}
//______________________________________________________________________________
{
// Return the length in centimeters of the current step
// TRACKR.ctrack = total curved path
- FlukaCallerCode_t caller = GetCaller();
- if (caller == kBXEntering || caller == kBXExiting ||
- caller == kENDRAW || caller == kUSDRAW ||
- caller == kUSTCKV || caller == kMGResumedTrack)
- return 0.0;
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kBXEntering || caller == kBXExiting ||
+ caller == kENDRAW || caller == kUSDRAW ||
+ caller == kUSTCKV || caller == kMGResumedTrack ||
+ caller == kSODRAW)
+ return 0.0;
else if (caller == kMGDRAW)
return TRACKR.ctrack;
else {
return TRACKR.cmtrck;
else if (caller == kMGResumedTrack)
return TRACKR.spausr[8];
+ else if (caller == kSODRAW)
+ return 0.0;
else {
- Warning("TrackLength", "track length not available");
+ Warning("TrackLength", "track length not available for caller %5d \n", caller);
return 0.0;
}
}
return TRACKR.atrack;
else if (caller == kMGResumedTrack)
return TRACKR.spausr[3];
+ else if (caller == kSODRAW) {
+ return (FLKSTK.agestk[FLKSTK.npflka]);
+ }
else {
Warning("TrackTime", "track time not available");
return 0.0;
// If coming from usdraw we just signal particle production - no edep
// If just first time after resuming, no edep for the primary
FlukaCallerCode_t caller = GetCaller();
+
if (caller == kBXExiting || caller == kBXEntering ||
- caller == kUSDRAW || caller == kMGResumedTrack) return 0.0;
+ caller == kUSDRAW || caller == kMGResumedTrack ||
+ caller == kSODRAW)
+ return 0.0;
Double_t sum = 0;
- if (TRACKR.mtrack > 1) printf("Edep: %6d\n", TRACKR.mtrack);
+ Int_t i = -1;
- for ( Int_t j=0;j<TRACKR.mtrack;j++) {
- sum +=TRACKR.dtrack[j];
- }
- if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
- return fRull + sum;
- else {
+ // Material with primary ionisation activated but number of primary electrons nprim = 0
+ if (fPrimaryElectronIndex == -2) return 0.0;
+ // nprim > 0
+ if ((i = fPrimaryElectronIndex) > -1) {
+ // Primary ionisation
+ sum = GetPrimaryElectronKineticEnergy(i);
+ if (sum > 100.) {
+ printf("edep > 100. %d %d %f \n", i, ALLDLT.nalldl, sum);
+ }
return sum;
+ } else {
+ // Normal ionisation
+ if (TRACKR.mtrack > 1) printf("Edep: %6d\n", TRACKR.mtrack);
+
+ for ( Int_t j=0;j<TRACKR.mtrack;j++) {
+ sum +=TRACKR.dtrack[j];
+ }
+ if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
+ return fRull + sum;
+ else {
+ return sum;
+ }
}
}
// and there is a call to endraw for energy deposition for each of them
// and they have the track number of their parent, but different identity (pdg)
// so we want to assign also their parent identity.
- if( (IsTrackStop() )
+
+ if( (IsTrackStop())
&& TRACKR.ispusr[mkbmx2 - 4] == TRACKR.ispusr[mkbmx2 - 1]
&& TRACKR.jtrack != TRACKR.ispusr[mkbmx2 - 3] ) {
if (fVerbosityLevel >=3)
<< " assign parent PDG=" << PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
return TRACKR.ispusr[mkbmx2 - 3]; // assign parent identity
}
- return TRACKR.jtrack;
+ if (TRACKR.jtrack <= 64){
+ return TRACKR.jtrack;
+ } else {
+ return TRACKR.j0trck;
+ }
}
// Return the id of the particle transported
// TRACKR.jtrack = identity number of the particle
FlukaCallerCode_t caller = GetCaller();
- if (caller != kEEDRAW) {
+ if (caller != kEEDRAW && caller != kSODRAW) {
return PDGFromId( CorrectFlukaId() );
}
+ else if (caller == kSODRAW) {
+ return PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
+ }
else
return -1000;
}
// Return charge of the track currently transported
// PAPROP.ichrge = electric charge of the particle
// TRACKR.jtrack = identity number of the particle
+
FlukaCallerCode_t caller = GetCaller();
- if (caller != kEEDRAW)
- return PAPROP.ichrge[CorrectFlukaId()+6];
+ if (caller != kEEDRAW && caller != kSODRAW)
+ return PAPROP.ichrge[CorrectFlukaId() + 6];
+ else if (caller == kSODRAW) {
+ Int_t ifl = PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
+ return PAPROP.ichrge[ifl + 6];
+ }
else
return -1000.0;
}
// PAPROP.am = particle mass in GeV
// TRACKR.jtrack = identity number of the particle
FlukaCallerCode_t caller = GetCaller();
- if (caller != kEEDRAW)
+ if (caller != kEEDRAW && caller != kSODRAW)
return PAPROP.am[CorrectFlukaId()+6];
+ else if (caller == kSODRAW) {
+ Int_t ifl = FLKSTK.iloflk[FLKSTK.npflka];
+ return PAPROP.am[ifl + 6];
+ }
else
return -1000.0;
}
Double_t TFluka::Etot() const
{
// TRACKR.etrack = total energy of the particle
- FlukaCallerCode_t caller = GetCaller();
- if (caller != kEEDRAW)
- return TRACKR.etrack;
- else
- return -1000.0;
+ FlukaCallerCode_t caller = GetCaller();
+ FlukaProcessCode_t icode = GetIcode();
+ if (caller != kEEDRAW && caller != kSODRAW && caller != kUSDRAW)
+ {
+ return TRACKR.etrack;
+ } else if (caller == kUSDRAW) {
+ if (icode == kEMFSCObrems ||
+ icode == kEMFSCOmoller ||
+ icode == kEMFSCObhabha ||
+ icode == kEMFSCOcompton ) {
+ return fPint[3];
+ }
+ else if (icode == kKASKADdray ||
+ icode == kKASKADbrems ||
+ icode == kKASKADpair) {
+ return (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);
+ }
+ }
+ else if (caller == kSODRAW) {
+ Int_t ist = FLKSTK.npflka;
+ Double_t p = FLKSTK.pmoflk[ist];
+ Int_t ifl = FLKSTK.iloflk[ist];
+ Double_t m = PAPROP.am[ifl + 6];
+ Double_t e = TMath::Sqrt(p * p + m * m);
+ return e;
+ }
+
+ return -1000.0;
}
//
//
// Return processes active in the current step
//
- FlukaProcessCode_t icode = GetIcode();
+ FlukaProcessCode_t icode = GetIcode();
+ FlukaCallerCode_t caller = GetCaller();
+
proc.Set(1);
TMCProcess iproc;
- switch (icode) {
- case kKASKADtimekill:
- case kEMFSCOtimekill:
- case kKASNEUtimekill:
- case kKASHEAtimekill:
- case kKASOPHtimekill:
- iproc = kPTOFlimit;
- break;
- case kKASKADstopping:
- case kKASKADescape:
- case kEMFSCOstopping1:
- case kEMFSCOstopping2:
- case kEMFSCOescape:
- case kKASNEUstopping:
- case kKASNEUescape:
- case kKASHEAescape:
- case kKASOPHescape:
- iproc = kPStop;
- break;
- case kKASOPHabsorption:
- iproc = kPLightAbsorption;
- break;
- case kKASOPHrefraction:
- iproc = kPLightRefraction;
- case kEMSCOlocaledep :
- iproc = kPPhotoelectric;
- break;
- default:
- iproc = ProdProcess(0);
+ if (caller == kBXEntering || caller == kBXExiting || caller == kEEDRAW || caller == kSODRAW) {
+ iproc = kPTransportation;
}
+ else if (caller == kUSTCKV) {
+ iproc = kPCerenkov;
+ } else {
+ switch (icode) {
+ case kEMFSCO:
+ if (Edep() > 0.) {
+ iproc = kPEnergyLoss;
+ } else {
+ iproc = kPTransportation;
+ }
+ break;
+ case kKASKAD:
+ if (Edep() > 0.) {
+ iproc = kPEnergyLoss;
+ } else {
+ iproc = kPTransportation;
+ }
+ break;
+ case kKASHEA:
+ case kKASNEU:
+ case kKASOPH:
+ case kKASKADescape:
+ case kEMFSCOescape:
+ case kKASNEUescape:
+ case kKASHEAescape:
+ case kKASOPHescape:
+ iproc = kPTransportation;
+ break;
+ case kKASKADtimekill:
+ case kEMFSCOtimekill:
+ case kKASNEUtimekill:
+ case kKASHEAtimekill:
+ case kKASOPHtimekill:
+ iproc = kPTOFlimit;
+ break;
+ case kKASKADstopping:
+ case kEMFSCOstopping1:
+ case kEMFSCOstopping2:
+ case kKASNEUstopping:
+ iproc = kPStop;
+ break;
+ case kKASOPHabsorption:
+ iproc = kPLightAbsorption;
+ break;
+ case kKASOPHrefraction:
+ iproc = kPLightRefraction;
+ break;
+ case kEMFSCOlocaldep :
+ iproc = kPPhotoelectric;
+ break;
+ default:
+ iproc = ProdProcess(0);
+ }
+ }
+
proc[0] = iproc;
return 1;
}
return fMCGeo->VolName(id);
}
+Int_t TFluka::MediumId(const Text_t* mediumName) const
+{
+ //
+ // Return the unique medium id for medium with name mediumName
+ TList *medlist = gGeoManager->GetListOfMedia();
+ TGeoMedium* med = (TGeoMedium*) medlist->FindObject(mediumName);
+ if (med) {
+ return (med->GetId());
+ } else {
+ return (-1);
+ }
+}
+
//______________________________________________________________________________
Int_t TFluka::VolId(const Text_t* volName) const
{
//
// Return the current volume name
//
- if (gGeoManager->IsOutside()) return 0;
+ if (gGeoManager->IsOutside()) return "OutOfWorld";
return gGeoManager->GetCurrentVolume()->GetName();
}
TDatabasePDG *pdgDB = TDatabasePDG::Instance();
- const Int_t kion=10000000;
-
const Double_t kAu2Gev = 0.9314943228;
const Double_t khSlash = 1.0545726663e-27;
const Double_t kErg2Gev = 1/1.6021773349e-3;
//
// Ions
//
-
pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,
- 0,3,"Ion",kion+10020);
+ 0,3,"Ion",GetIonPdg(1,2));
pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,
- khShGev/(12.33*kYear2Sec),3,"Ion",kion+10030);
+ khShGev/(12.33*kYear2Sec),3,"Ion",GetIonPdg(1,3));
pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,
- khShGev/(12.33*kYear2Sec),6,"Ion",kion+20040);
+ khShGev/(12.33*kYear2Sec),6,"Ion",GetIonPdg(2,4));
pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE,
- 0,6,"Ion",kion+20030);
+ 0,6,"Ion",GetIonPdg(2,3));
}
//
}
//______________________________________________________________________________
-Double_t TFluka::GetPrimaryElectronKineticEnergy(Int_t i)
+Double_t TFluka::GetPrimaryElectronKineticEnergy(Int_t i) const
{
- Double_t ekin = -1.;
// Returns kinetic energy of primary electron i
+
+ Double_t ekin = -1.;
+
if (i >= 0 && i < ALLDLT.nalldl) {
ekin = ALLDLT.talldl[i];
} else {
}
return ekin;
}
+
+void TFluka::GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Double_t& z) const
+{
+ // Returns position of primary electron i
+ if (i >= 0 && i < ALLDLT.nalldl) {
+ x = ALLDLT.xalldl[i];
+ y = ALLDLT.yalldl[i];
+ z = ALLDLT.zalldl[i];
+ return;
+ } else {
+ Warning("GetPrimaryElectronPosition",
+ "Primary electron index out of range %d %d \n",
+ i, ALLDLT.nalldl);
+ return;
+ }
+ return;
+}
+
+Int_t TFluka::GetIonPdg(Int_t z, Int_t a, Int_t i) const
+{
+// Acording to
+// http://cepa.fnal.gov/psm/stdhep/pdg/montecarlorpp-2006.pdf
+
+ return 1000000000 + 10*1000*z + 10*a + i;
+}
+
+void TFluka::PrimaryIonisationStepping(Int_t nprim)
+{
+// Call Stepping for primary ionisation electrons
+ Int_t i;
+// Protection against nprim > mxalld
+
+// Multiple steps for nprim > 0
+ if (nprim > 0) {
+ for (i = 0; i < nprim; i++) {
+ SetCurrentPrimaryElectronIndex(i);
+ (TVirtualMCApplication::Instance())->Stepping();
+ if (i == 0) SetTrackIsNew(kFALSE);
+ }
+ } else {
+ // No primary electron ionisation
+ // Call Stepping anyway but flag nprim = 0 as index = -2
+ SetCurrentPrimaryElectronIndex(-2);
+ (TVirtualMCApplication::Instance())->Stepping();
+ }
+ // Reset the index
+ SetCurrentPrimaryElectronIndex(-1);
+}
+
+//______________________________________________________________________
+Float_t* TFluka::CreateFloatArray(Double_t* array, Int_t size) const
+{
+// Converts Double_t* array to Float_t*,
+// !! The new array has to be deleted by user.
+// ---
+
+ Float_t* floatArray;
+ if (size>0) {
+ floatArray = new Float_t[size];
+ for (Int_t i=0; i<size; i++)
+ if (array[i] >= FLT_MAX )
+ floatArray[i] = FLT_MAX/100.;
+ else
+ floatArray[i] = array[i];
+ }
+ else {
+ //floatArray = 0;
+ floatArray = new Float_t[1];
+ }
+ return floatArray;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff