#include <TList.h>
#include "TFluka.h"
+#include "TFlukaIon.h"
#include "TFlukaCodes.h"
#include "TCallf77.h" //For the fortran calls
#include "Fdblprc.h" //(DBLPRC) fluka common
fXsco(0),
fYsco(0),
fZsco(0),
+ fPItime(0),
+ fPIlength(0),
+ fNPI(0),
fTrackIsEntering(kFALSE),
fTrackIsExiting(kFALSE),
fTrackIsNew(kFALSE),
fStopEvent(kFALSE),
fStopRun(kFALSE),
fPrimaryElectronIndex(-1),
+ fLowEnergyNeutronTransport(kFALSE),
fMaterials(0),
fNVolumes(0),
fCurrentFlukaRegion(-1),
fGeom(0),
fMCGeo(0),
fUserConfig(0),
- fUserScore(0)
+ fUserScore(0),
+ fUserIon(0)
{
//
// Default constructor
fXsco(0),
fYsco(0),
fZsco(0),
+ fPItime(0),
+ fPIlength(0),
+ fNPI(0),
fTrackIsEntering(kFALSE),
fTrackIsExiting(kFALSE),
fTrackIsNew(kFALSE),
fStopEvent(kFALSE),
fStopRun(kFALSE),
fPrimaryElectronIndex(-1),
+ fLowEnergyNeutronTransport(kFALSE),
fMaterials(0),
fNVolumes(0),
fCurrentFlukaRegion(-1),
fGeom(0),
fMCGeo(0),
fUserConfig(new TObjArray(100)),
- fUserScore(new TObjArray(100))
+ fUserScore(new TObjArray(100)),
+ fUserIon(0)
{
// create geometry interface
for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
cout << "<== TFluka::~TFluka() destructor called." << endl;
if (fMaterials) delete [] fMaterials;
- delete fGeom;
- delete fMCGeo;
+// delete fGeom;
+// delete fMCGeo;
if (fUserConfig) {
fUserConfig->Delete();
// Add ions to PDG Data base
//
AddParticlesToPdgDataBase();
- //
+ //
}
// 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;
}
- printf("Mixture (2) %5d %10s %5d \n", kmat, name, nlmat);
gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
}
strncmp(param, "HADR", 4) == 0 ||
strncmp(param, "LOSS", 4) == 0 ||
strncmp(param, "MULS", 4) == 0 ||
- strncmp(param, "RAYL", 4) == 0)
+ strncmp(param, "RAYL", 4) == 0 ||
+ strncmp(param, "STRA", 4) == 0)
{
process = kTRUE;
}
//
//
// 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
+//
+// npckov: number of sampling points
+// ppckov: energy values
+// absco: absorption length
+// effic: quantum efficiency
+// rindex: refraction index
//
-// Double_t version not implemented
+
+//
+// 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) {
//
-// // Double_t version not implemented
+// Set Cerenkov properties for medium itmed
+//
+// 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 !!!
+ }
}
//____________________________________________________________________________
//_____________________________________________________________________________
Int_t TFluka::IdFromPDG(Int_t pdg) const
{
+
//
// Return Fluka code from PDG and pseudo ENDF code
-
+ Int_t idSpecial[4] = {TFlukaIon::GetIonPdg(2,4),
+ TFlukaIon::GetIonPdg(2,3),
+ TFlukaIon::GetIonPdg(1,3),
+ TFlukaIon::GetIonPdg(1,2)};
// Catch the feedback photons
if (pdg == 50000051) return (kFLUKAoptical);
+ // Ion as primary
+ for (Int_t i = 0; i < 4; i++) {
+ if (pdg == idSpecial[i]) return (i + kFLUKAcodemin);
+ }
+
+ if ((!fUserIon && pdg == TFlukaIon::GetIonPdg(6,12)) ||
+ ( fUserIon && pdg == fUserIon->GetPdgCode()))
+ return (-2);
+
// MCIHAD() goes from pdg to fluka internal.
Int_t intfluka = mcihad(pdg);
// KPTOIP array goes from internal to official
{
//
// Return PDG code and pseudo ENDF code from Fluka code
- // Alpha He3 Triton Deuteron gen. ion opt. photon
- Int_t idSpecial[6] = {GetIonPdg(2,4), GetIonPdg(2, 3), GetIonPdg(1,3), GetIonPdg(1,2), GetIonPdg(0,0), 50000050};
+ Int_t idSpecial[6] = {TFlukaIon::GetIonPdg(2,4), // alpha
+ TFlukaIon::GetIonPdg(2,3), // He3
+ TFlukaIon::GetIonPdg(1,3), // triton
+ TFlukaIon::GetIonPdg(1,2), // deuteron
+ TFlukaIon::GetIonPdg(0,0), // gen. ion
+ 50000050};
// IPTOKP array goes from official to internal
if (id == kFLUKAoptical) {
//
// Construct file names
FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
- TString sFlukaVmcCoreInp = getenv("ALICE_ROOT");
- sFlukaVmcCoreInp +="/TFluka/input/";
TString sFlukaVmcTmp = "flukaMat.inp";
TString sFlukaVmcInp = GetInputFileName();
- sFlukaVmcCoreInp += GetCoreInputFileName();
+ TString sFlukaVmcCoreInp = GetCoreInputFileName();
// Open files
if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
// Add RANDOMIZ card
fprintf(pFlukaVmcInp,"RANDOMIZ %10.1f%10.0f\n", 1., Float_t(gRandom->GetSeed()));
+// User defined ion
+ if (fUserIon) fUserIon->WriteUserInputCard(pFlukaVmcInp);
// Add START and STOP card
fprintf(pFlukaVmcInp,"START %10.1f\n",fEventsPerRun);
fprintf(pFlukaVmcInp,"STOP \n");
// Int_t mreg=0, latt=0;
// fGeom->GetCurrentRegion(mreg, latt);
+
+
Int_t mreg = fGeom->GetCurrentRegion();
STEPSZ.stepmx[mreg - 1] = step;
}
Int_t i = -1;
if ((i = fPrimaryElectronIndex) > -1) {
// Primary Electron Ionisation
- Double_t x, y, z;
- GetPrimaryElectronPosition(i, x, y, z);
+ Double_t x, y, z, t;
+ GetPrimaryElectronPosition(i, x, y, z, t);
position.SetX(x);
position.SetY(y);
position.SetZ(z);
- position.SetT(TRACKR.atrack);
+ position.SetT(t);
} else {
position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
else if (caller == kMGDRAW) {
Int_t i = -1;
if ((i = fPrimaryElectronIndex) > -1) {
- GetPrimaryElectronPosition(i, x, y, z);
+ Double_t t;
+ GetPrimaryElectronPosition(i, x, y, z, t);
} else {
x = TRACKR.xtrack[TRACKR.ntrack];
y = TRACKR.ytrack[TRACKR.ntrack];
momentum.SetPz(p * FLKSTK.tzflk[ist]);
momentum.SetE(e);
} else if (caller == kUSDRAW) {
- if (icode == 208 || icode == 210 || icode == 212 || icode == 219) {
+ 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));
+ 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);
py = p * FLKSTK.tyflk[ist];
pz = p * FLKSTK.tzflk[ist];
} else if (caller == kUSDRAW) {
- if (icode == 208 || icode == 210 || icode == 212 || icode == 219) {
+ 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;
// 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 ||
- caller == kSODRAW)
+ if (caller == kMGDRAW) {
+ Int_t i;
+ if ((i = fPrimaryElectronIndex) > -1) {
+ if (i > 0) {
+ return (fPIlength[i] - fPIlength[i-1]);
+ } else {
+ Double_t s (TRACKR.ctrack - (fPIlength[fNPI - 1] - fPIlength[0]));
+ return s;
+ }
+ } else {
+ return TRACKR.ctrack;
+ }
+ } else 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 {
- Warning("TrackStep", "track step not available");
- return 0.0;
- }
+ } else {
+ Warning("TrackStep", "track step not available");
+ return 0.0;
+ }
}
//______________________________________________________________________________
{
// TRACKR.cmtrck = cumulative curved path since particle birth
FlukaCallerCode_t caller = GetCaller();
- if (caller == kBXEntering || caller == kBXExiting ||
- caller == kENDRAW || caller == kUSDRAW || caller == kMGDRAW ||
- caller == kUSTCKV)
- return TRACKR.cmtrck;
+ if (caller == kMGDRAW) {
+ Int_t i;
+ if ((i = fPrimaryElectronIndex) > -1) {
+ return fPIlength[i];
+ } else {
+ return TRACKR.cmtrck;
+ }
+
+ } else if (caller == kBXEntering || caller == kBXExiting ||
+ caller == kENDRAW || caller == kUSDRAW || caller == kUSTCKV)
+ return TRACKR.cmtrck;
else if (caller == kMGResumedTrack)
- return TRACKR.spausr[8];
+ return TRACKR.spausr[8];
else if (caller == kSODRAW)
return 0.0;
else {
- Warning("TrackLength", "track length not available for caller %5d \n", caller);
- return 0.0;
- }
+ Warning("TrackLength", "track length not available for caller %5d \n", caller);
+ return 0.0;
+ }
}
+
//______________________________________________________________________________
Double_t TFluka::TrackTime() const
{
// Return the current time of flight of the track being transported
// TRACKR.atrack = age of the particle
FlukaCallerCode_t caller = GetCaller();
- if (caller == kBXEntering || caller == kBXExiting ||
- caller == kENDRAW || caller == kUSDRAW || caller == kMGDRAW ||
- caller == kUSTCKV)
+ if (caller == kMGDRAW) {
+ Int_t i;
+ if ((i = fPrimaryElectronIndex) > -1) {
+ Double_t t = fPItime[i];
+ return t;
+ } else {
+ return TRACKR.atrack;
+ }
+ } else if (caller == kBXEntering || caller == kBXExiting ||
+ caller == kENDRAW || caller == kUSDRAW ||
+ caller == kUSTCKV)
return TRACKR.atrack;
else if (caller == kMGResumedTrack)
return TRACKR.spausr[3];
Double_t TFluka::Etot() const
{
// TRACKR.etrack = total energy of the particle
- FlukaCallerCode_t caller = GetCaller();
- if (caller != kEEDRAW && caller != kSODRAW)
- return TRACKR.etrack;
+ 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 {
+ return TRACKR.etrack;
+ }
+
+ }
else if (caller == kSODRAW) {
Int_t ist = FLKSTK.npflka;
Double_t p = FLKSTK.pmoflk[ist];
Double_t e = TMath::Sqrt(p * p + m * m);
return e;
}
- else
- return -1000.0;
+ printf("Etot %5d %5d \n", caller, icode);
+
+ return -1000.0;
}
//
//______________________________________________________________________________
Bool_t TFluka::IsTrackAlive() const
{
-// means not disappeared or not out
- if (IsTrackDisappeared() || IsTrackOut() ) return 0;
- else return 1;
+// Means not disappeared or not out
+ FlukaProcessCode_t icode = GetIcode();
+
+ if (IsTrackOut() ||
+ IsTrackStop() ||
+ icode == kKASKADinelint || // inelastic interaction
+ icode == kKASKADdecay || // particle decay
+ icode == kEMFSCOanniflight || // in-flight annihilation
+ icode == kEMFSCOannirest || // annihilation at rest
+ icode == kEMFSCOpair || // pair production
+ icode == kEMFSCOphotoel || // Photoelectric effect
+ icode == kKASNEUhadronic // hadronic interaction
+ )
+ {
+ // Exclude the cases for which the particle has disappeared (paused) but will reappear later (= alive).
+ return 0;
+ } else {
+ return 1;
+ }
}
//
//
FlukaProcessCode_t icode = GetIcode();
FlukaCallerCode_t caller = GetCaller();
-
proc.Set(1);
TMCProcess iproc;
- if (caller == kBXEntering || caller == kBXExiting || caller == kEEDRAW) {
+ if (caller == kBXEntering || caller == kBXExiting || caller == kEEDRAW || caller == kSODRAW) {
iproc = kPTransportation;
+ }
+ else if (caller == kUSTCKV) {
+ iproc = kPCerenkov;
} else {
switch (icode) {
case kEMFSCO:
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:
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 kKASKADinelint:
+ case kKASNEUhadronic:
+ iproc = kPHadronic;
+ break;
+ case kKASKADinelarecoil:
+ iproc = kPHadronic;
+ break;
+ case kKASKADnelint:
+ iproc = kPHElastic;
break;
case kKASOPHabsorption:
iproc = kPLightAbsorption;
break;
case kKASOPHrefraction:
iproc = kPLightRefraction;
+ break;
case kEMFSCOlocaldep :
iproc = kPPhotoelectric;
break;
//
// Return the current volume name
//
- if (gGeoManager->IsOutside()) return "Outside FLUKA Geometry !";
+ if (gGeoManager->IsOutside()) return "OutOfWorld";
return gGeoManager->GetCurrentVolume()->GetName();
}
// Ions
//
pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,
- 0,3,"Ion",GetIonPdg(1,2));
+ 0,3,"Ion",TFlukaIon::GetIonPdg(1,2));
pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,
- khShGev/(12.33*kYear2Sec),3,"Ion",GetIonPdg(1,3));
+ khShGev/(12.33*kYear2Sec),3,"Ion",TFlukaIon::GetIonPdg(1,3));
pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,
- khShGev/(12.33*kYear2Sec),6,"Ion",GetIonPdg(2,4));
+ khShGev/(12.33*kYear2Sec),6,"Ion",TFlukaIon::GetIonPdg(2,4));
pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE,
- 0,6,"Ion",GetIonPdg(2,3));
+ 0,6,"Ion",TFlukaIon::GetIonPdg(2,3));
+//
+// Default user ion
+ TFlukaIon::AddIon(12, 6);
+
+//
+//
+//
+// Special particles
+//
+ pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,
+ 0,0,"Special",GetSpecialPdg(50));
+ pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
+ 0,0,"Special",GetSpecialPdg(51));
}
+
//
// Info about primary ionization electrons
//
return ekin;
}
-void TFluka::GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Double_t& z) const
+void TFluka::GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Double_t& z, Double_t& t) 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];
+ t = ALLDLT.talldl[i];
return;
} else {
Warning("GetPrimaryElectronPosition",
return;
}
-Int_t TFluka::GetIonPdg(Int_t z, Int_t a, Int_t i) const
+
+//__________________________________________________________________
+Int_t TFluka::GetSpecialPdg(Int_t number) const
{
-// Acording to
-// http://cepa.fnal.gov/psm/stdhep/pdg/montecarlorpp-2006.pdf
+// Numbering for special particles
+
+ return 50000000 + number;
+}
- 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
+ Int_t i;
+ fNPI = nprim;
if (nprim > 0) {
+ CalcPrimaryIonisationTime();
for (i = 0; i < nprim; i++) {
SetCurrentPrimaryElectronIndex(i);
(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;
+}
+
+void TFluka::CalcPrimaryIonisationTime()
+{
+ // Calculates the primary ionisation time
+ if (fPItime) delete [] fPItime;
+ fPItime = new Double_t[fNPI];
+ if (fPIlength) delete [] fPIlength;
+ fPIlength = new Double_t[fNPI];
+ //
+ Double_t px, py, pz, e, t;
+ TrackMomentum(px, py, pz, e);
+ Double_t p = TMath::Sqrt(px * px + py * py + pz * pz);
+ Double_t beta = p / e;
+ Double_t x0, y0, z0;
+ fPItime[fNPI -1] = TRACKR.atrack;
+ fPIlength[fNPI -1] = TRACKR.cmtrck;
+ GetPrimaryElectronPosition(fNPI - 1, x0, y0, z0, t);
+ if (fNPI > 1) {
+ for (Int_t i = fNPI - 2; i > -1; i--) {
+ Double_t x, y, z, t;
+ GetPrimaryElectronPosition(i, x, y, z, t);
+ Double_t ds = TMath::Sqrt((x-x0) * (x-x0) + (y-y0) * (y-y0) + (z-z0) * (z-z0));
+ fPItime[i] = fPItime[i+1] - ds / (beta * 2.99792458e10);
+ fPIlength[i] = fPIlength[i+1] - ds;
+ x0 = x; y0 = y; z0 = z;
+ }
+ }
+
+}
+
+Bool_t TFluka::DefineIon(const char* name , Int_t z, Int_t a, Int_t q, Double_t exE, Double_t mass)
+{
+ // User defined ion that can be used as a primary
+ if (fUserIon) {
+ Warning("DefineIon", "Only one user ion can be defined !");
+ return kFALSE;
+ } else {
+ fUserIon = new TFlukaIon(name, z, a, q, exE, mass);
+ return kTRUE;
+ }
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff