- Improved IsAlive definition.
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 e = TMath::Sqrt(p * p + m * m);
return e;
}
+ printf("Etot %5d %5d \n", caller, icode);
return -1000.0;
}
//
FlukaProcessCode_t icode = GetIcode();
FlukaCallerCode_t caller = GetCaller();
-
proc.Set(1);
TMCProcess iproc;
if (caller == kBXEntering || caller == kBXExiting || caller == kEEDRAW || caller == kSODRAW) {
case kEMFSCOstopping2:
case kKASNEUstopping:
iproc = kPStop;
+ break;
+ case kKASKADinelint:
+ case kKASNEUhadronic:
+ iproc = kPHadronic;
+ break;
+ case kKASKADinelarecoil:
+ iproc = kPHadronic;
+ break;
+ case kKASKADnelint:
+ iproc = kPHElastic;
break;
case kKASOPHabsorption:
iproc = kPLightAbsorption;
0,0,"Special",GetSpecialPdg(51));
}
+void TFluka::AddIon(Int_t a, Int_t z) const
+{
+
+ // Add a new ion
+ TDatabasePDG *pdgDB = TDatabasePDG::Instance();
+ const Double_t kAu2Gev = 0.9314943228;
+ Int_t pdg = GetIonPdg(z, a);
+ if (pdgDB->GetParticle(pdg)) return;
+
+ pdgDB->AddParticle(Form("Iion A = %5d Z = %5d", a, z),"Ion", Float_t(a) * kAu2Gev + 8.071e-3, kTRUE,
+ 0, 3 * z, "Ion", pdg);
+}
+
//
// Info about primary ionization electrons
//
void GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Double_t& z, Double_t& t) const;
void SetCurrentPrimaryElectronIndex(Int_t i) {fPrimaryElectronIndex = i;}
void PrimaryIonisationStepping(Int_t nprim);
+
+ void AddIon(Int_t a, Int_t z) const;
+ Int_t GetIonPdg(Int_t z, Int_t a, Int_t i = 0) const;
private:
// Copy constructor and operator= declared but not implemented (-Weff++ flag)
void PrintHeader();
void AddParticlesToPdgDataBase() const;
- Int_t GetIonPdg(Int_t z, Int_t a, Int_t i = 0) const;
+
Int_t GetSpecialPdg(Int_t number) const;
Float_t* CreateFloatArray(Double_t* array, Int_t size) const;
}
Int_t verbosityLevel = fluka->GetVerbosityLevel();
-
if (TRACKR.jtrack < -6) {
// from -7 to -12 = "heavy" fragment
// assing parent id
cout << "mgdraw: (heavy fragment) jtrack < -6 =" << TRACKR.jtrack
<< " assign parent pdg=" << fluka->PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
}
- TRACKR.jtrack = TRACKR.ispusr[mkbmx2 - 3];
+// TRACKR.jtrack = TRACKR.ispusr[mkbmx2 - 3];
}
cppstack->SetCurrentTrack(trackId);
#include "Fflkstk.h" //(FLKSTK) fluka common
#include "Ftrackr.h" //(TRACKR) fluka common
#include "Fgenstk.h" //(GENSTK) fluka common
+#include "Ffheavy.h" //(FHEAVY) fluka common
//Virtual MC
//* SeT User PRoperties for Fluka particles *
//* *
//*----------------------------------------------------------------------*
+// Get the pointer to the VMC
+ TFluka* fluka = (TFluka*) gMC;
// FLKSTK.npflka = stack pointer
// FLKSTK.louse = user flag
FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 4] = TRACKR.ispusr[mkbmx2 - 1]; // current track number
FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 5] = npprmr; // flag npprmr>0
-// Get the pointer to the VMC
- TFluka* fluka = (TFluka*) gMC;
+
Int_t verbosityLevel = fluka->GetVerbosityLevel();
Bool_t debug = (verbosityLevel>=3)? kTRUE:kFALSE;
Int_t parent = TRACKR.ispusr[mkbmx2-1];
Int_t kpart = GENSTK.kpart[numsec-1];
- if (kpart < -6) return; // -7 to -12 = "heavy" fragment
- Int_t pdg = fluka->PDGFromId(kpart);
-
+ Int_t pdg;
+ Int_t a, z;
+
+ if (kpart < -6) {
+ kpart = -kpart;
+ z = kpart / 100000;
+ a = kpart - z * 100000;
+ a /= 100;
+ pdg = fluka->GetIonPdg(z, a);
+ fluka->AddIon(a, z);
+ } else {
+ pdg = fluka->PDGFromId(kpart);
+ }
+
Double_t px = GENSTK.plr[numsec-1] * GENSTK.cxr[numsec-1];
Double_t py = GENSTK.plr[numsec-1] * GENSTK.cyr[numsec-1];
Double_t pz = GENSTK.plr[numsec-1] * GENSTK.czr[numsec-1];
- Double_t e = GENSTK.tki[numsec-1] + PAPROP.am[GENSTK.kpart[numsec-1]+6];
-
+ Double_t e;
+
+ if (kpart < 100000) {
+ e = GENSTK.tki[numsec-1] + PAPROP.am[GENSTK.kpart[numsec-1]+6];
+ } else {
+ e = GENSTK.tki[numsec-1] + Double_t(a) * 0.93149 + 8.071e-3;
+ }
+
+
Double_t vx = xx;
Double_t vy = yy;
Double_t vz = zz;
Double_t poly = GENSTK.cyrpol[numsec-1];
Double_t polz = GENSTK.czrpol[numsec-1];
-
-
TMCProcess mech = kPHadronic;
if (EVTFLG.ldecay == 1) {
mech = kPDecay;
vx, vy, vz, tof,
polx, poly, polz,
mech, ntr, weight, is);
- if (debug)
+ if (debug) {
+
cout << endl << " !!! stuprf: ntr=" << ntr << " pdg " << pdg << " parent=" << parent
<< " parent_pdg="<< fluka->PDGFromId(TRACKR.jtrack) << " numsec "
<< numsec << " npprmr " << npprmr << " icode=" << fluka->GetIcode()
<< "kin. energy [keV] = " << GENSTK.tki[numsec-1] * 1.e6 << endl << endl;
+ }
+
//
// Save current track number