//
#include <Riostream.h>
+#include <TList.h>
#include "TFluka.h"
#include "TFlukaCodes.h"
fStopped(kFALSE),
fStopEvent(kFALSE),
fStopRun(kFALSE),
+ fPrimaryElectronIndex(-1),
fMaterials(0),
fNVolumes(0),
fCurrentFlukaRegion(-1),
fStopped(kFALSE),
fStopEvent(kFALSE),
fStopRun(kFALSE),
+ fPrimaryElectronIndex(-1),
fMaterials(0),
fNVolumes(0),
fCurrentFlukaRegion(-1),
}
fApplication->InitGeometry();
-
+ fApplication->ConstructOpGeometry();
//
// Add ions to PDG Data base
//
AddParticlesToPdgDataBase();
+ //
+
+
}
// 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);
}
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]);
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];
+ 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 == kMGResumedTrack) {
x = TRACKR.spausr[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;
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;
+ }
}
break;
case kKASOPHrefraction:
iproc = kPLightRefraction;
- case kEMSCOlocaledep :
+ case kEMFSCOlocaldep :
iproc = kPPhotoelectric;
break;
default:
}
//______________________________________________________________________________
-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;
+}
+
+