X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TFluka%2FTFluka.cxx;h=0fd9649ab1648ca2fcad8b974d17d1f9a1442b43;hb=81842c1640a56ed6be5ead79dfb0245de1f45904;hp=cd7a3d0490d3bb7a2bfb560d6c1fd7903951003e;hpb=09cdde8a9fdb9b548b82e3a0598be337608ee0c1;p=u%2Fmrichter%2FAliRoot.git

diff --git a/TFluka/TFluka.cxx b/TFluka/TFluka.cxx
index cd7a3d0490d..0fd9649ab16 100644
--- a/TFluka/TFluka.cxx
+++ b/TFluka/TFluka.cxx
@@ -30,6 +30,7 @@
 //
 
 #include <Riostream.h>
+#include <TList.h>
 
 #include "TFluka.h"
 #include "TFlukaCodes.h"
@@ -137,6 +138,7 @@ TFluka::TFluka()
    fStopped(kFALSE),
    fStopEvent(kFALSE),
    fStopRun(kFALSE),
+   fPrimaryElectronIndex(-1),
    fMaterials(0),
    fNVolumes(0),
    fCurrentFlukaRegion(-1),
@@ -174,6 +176,7 @@ TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupporte
    fStopped(kFALSE),
    fStopEvent(kFALSE),
    fStopRun(kFALSE),
+   fPrimaryElectronIndex(-1),
    fMaterials(0),
    fNVolumes(0),
    fCurrentFlukaRegion(-1),
@@ -247,11 +250,14 @@ void TFluka::Init() {
     }
 
     fApplication->InitGeometry();
-
+    fApplication->ConstructOpGeometry();
     //
     // Add ions to PDG Data base
     //
      AddParticlesToPdgDataBase();
+     //
+
+     
 }
 
 
@@ -995,7 +1001,7 @@ Int_t TFluka::PDGFromId(Int_t id) const
 // 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    
@@ -1116,21 +1122,22 @@ Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
 }
 
 
-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);
 }
 
@@ -1335,11 +1342,22 @@ void TFluka::TrackPosition(TLorentzVector& position) const
     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]);
@@ -1369,14 +1387,19 @@ void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
   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];
@@ -1550,18 +1573,34 @@ Double_t TFluka::Edep() const
   // 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;
+      }
   }
 }
 
@@ -1572,7 +1611,8 @@ Int_t TFluka::CorrectFlukaId() const
    // 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)
@@ -1582,7 +1622,11 @@ Int_t TFluka::CorrectFlukaId() const
                << " 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;
+   }
 }
 
 
@@ -1900,7 +1944,7 @@ Int_t TFluka::StepProcesses(TArrayI &proc) const
         break;
     case kKASOPHrefraction:
         iproc = kPLightRefraction;
-    case kEMSCOlocaledep : 
+    case kEMFSCOlocaldep : 
         iproc = kPPhotoelectric;
         break;
     default:
@@ -2289,10 +2333,11 @@ Int_t TFluka::GetNPrimaryElectrons()
 }
 
 //______________________________________________________________________________
-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 {
@@ -2302,3 +2347,22 @@ Double_t TFluka::GetPrimaryElectronKineticEnergy(Int_t i)
     }
     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;
+}
+
+