Correct loop over primary ionisations in case of resuming transport.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 5 Oct 2007 20:00:27 +0000 (20:00 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 5 Oct 2007 20:00:27 +0000 (20:00 +0000)
TFluka/TFluka.cxx
TFluka/TFluka.h
TFluka/mgdraw.cxx

index 22e0d6b..cde7661 100644 (file)
@@ -655,8 +655,7 @@ void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
        strncmp(param, "HADR",  4) == 0 ||
        strncmp(param, "LOSS",  4) == 0 ||
        strncmp(param, "MULS",  4) == 0 ||
-       strncmp(param, "RAYL",  4) == 0 ||
-       strncmp(param, "STRA",  4) == 0) 
+       strncmp(param, "RAYL",  4) == 0) 
    {
        process = kTRUE;
    } 
@@ -1758,7 +1757,6 @@ Bool_t   TFluka::IsTrackDisappeared() const
 {
 // All inelastic interactions and decays
 // fIcode from usdraw
-    
   FlukaProcessCode_t icode = GetIcode();
   if (icode == kKASKADinelint    || // inelastic interaction
       icode == kKASKADdecay      || // particle decay
@@ -1801,27 +1799,8 @@ Bool_t   TFluka::IsTrackStop() const
 Bool_t   TFluka::IsTrackAlive() const
 {
 // means not disappeared or not out
-    FlukaProcessCode_t icode = GetIcode();
-
-    if (IsTrackOut() || IsTrackStop()) {
-       return 0;
-    } else if 
-       (
-       IsTrackDisappeared()    &&
-       icode != kKASKADdray    &&
-       icode != kKASKADpair    &&
-       icode != kKASKADbrems   &&
-       icode != kEMFSCObrems   &&
-       icode != kEMFSCOmoller  &&
-       icode != kEMFSCObhabha  &&
-       icode != kEMFSCOcompton
-       ) 
-    {
-       // Exclude the cases for which the particle has disappeared (paused) but will reappear later (= alive).
-       return 0;
-    } else {
-       return 1;
-    }
+  if (IsTrackDisappeared() || IsTrackOut() ) return 0;
+  else return 1;
 }
 
 //
@@ -2358,8 +2337,10 @@ Double_t TFluka::GetPrimaryElectronKineticEnergy(Int_t i) const
     // Returns kinetic energy of primary electron i
 
     Double_t ekin = -1.;
+    
     if (i >= 0 && i < ALLDLT.nalldl) {
-        ekin =  ALLDLT.talldl[i];
+       Int_t j = ALLDLT.nalldl - 1 - i;
+        ekin =  ALLDLT.talldl[j];
     } else {
         Warning("GetPrimaryElectronKineticEnergy",
                 "Primary electron index out of range %d %d \n",
@@ -2372,9 +2353,10 @@ void TFluka::GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Doubl
 {
     // Returns position  of primary electron i
         if (i >= 0 && i < ALLDLT.nalldl) {
-           x = ALLDLT.xalldl[i];
-           y = ALLDLT.yalldl[i];
-           z = ALLDLT.zalldl[i];
+           Int_t j = ALLDLT.nalldl - 1 - i;
+           x = ALLDLT.xalldl[j];
+           y = ALLDLT.yalldl[j];
+           z = ALLDLT.zalldl[j];
            return;
        } else {
            Warning("GetPrimaryElectronPosition",
@@ -2393,3 +2375,25 @@ Int_t TFluka::GetIonPdg(Int_t z, Int_t a, Int_t i) const
   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);
+}
index 314336f..3ef6f42 100644 (file)
@@ -400,7 +400,7 @@ class TFluka : public TVirtualMC {
   Double_t GetPrimaryElectronKineticEnergy(Int_t i) const;
   void     GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Double_t& z) const;
   void     SetCurrentPrimaryElectronIndex(Int_t i)  {fPrimaryElectronIndex = i;}
-
+  void     PrimaryIonisationStepping(Int_t nprim);  
  private:
    
   // Copy constructor and operator= declared but not implemented (-Weff++ flag)
@@ -410,9 +410,8 @@ class TFluka : public TVirtualMC {
   void  PrintHeader();
   void  AddParticlesToPdgDataBase() const;
   Int_t GetIonPdg(Int_t z, Int_t a, Int_t i = 0) const;
-  
-  //
 
+  //
   Int_t   fVerbosityLevel; //Verbosity level (0 lowest - 3 highest)
   Int_t   fNEvent;         //Current event number
   TString fInputFileName;     //Name of the real input file 
index 33accf8..6490139 100644 (file)
@@ -104,32 +104,16 @@ void mgdraw(Int_t& icode, Int_t& mreg)
 //
 // Primary ionsiations
            Int_t nprim = ALLDLT.nalldl;
-// Protection against nprim > mxalld
            if (nprim >= mxalld) {
                nprim = mxalld;
                Warning("mgdraw", "nprim > mxalld, nprim: %6d  pdg: %6d mreg %6d p %13.3f step %13.3f\n", 
-                      ALLDLT.nalldl, 
-                      fluka->PDGFromId(TRACKR.jtrack), 
-                      mreg, 
-                      TRACKR.ptrack, 
-                      TRACKR.ctrack);
-               
-           }
-// Multiple steps for nprim > 0
-           if (nprim > 0) {
-               for (Int_t i = 0; i < nprim; i++) {
-                   fluka->SetCurrentPrimaryElectronIndex(i);
-                   (TVirtualMCApplication::Instance())->Stepping();
-                   if (i == 0) fluka->SetTrackIsNew(kFALSE);
-               }
-           } else {
-               // No primary electron ionisation
-               // Call Stepping anyway but flag nprim = 0 as index = -2
-               fluka->SetCurrentPrimaryElectronIndex(-2);
-               (TVirtualMCApplication::Instance())->Stepping();
+                       ALLDLT.nalldl, 
+                       fluka->PDGFromId(TRACKR.jtrack), 
+                       mreg, 
+                       TRACKR.ptrack, 
+                       TRACKR.ctrack);
            }
-           // Reset the index
-           fluka->SetCurrentPrimaryElectronIndex(-1);
+           fluka->PrimaryIonisationStepping(nprim);
        } else {
            // Single step
            (TVirtualMCApplication::Instance())->Stepping();
@@ -160,10 +144,23 @@ void mgdraw(Int_t& icode, Int_t& mreg)
 
         fluka->SetTrackIsNew(kFALSE);
         fluka->SetCaller(kMGDRAW);
-       if (msd > 0) fluka->SetCurrentPrimaryElectronIndex(-2);
-        (TVirtualMCApplication::Instance())->Stepping();
-       if (msd > 0) fluka->SetCurrentPrimaryElectronIndex(-1);
-    }
+       if (msd == 0) {
+           (TVirtualMCApplication::Instance())->Stepping();
+       } else {
+           Int_t nprim = ALLDLT.nalldl;
+// Protection against nprim > mxalld
+           if (nprim >= mxalld) {
+               nprim = mxalld;
+               Warning("mgdraw", "nprim > mxalld, nprim: %6d  pdg: %6d mreg %6d p %13.3f step %13.3f\n", 
+                       ALLDLT.nalldl, 
+                       fluka->PDGFromId(TRACKR.jtrack), 
+                       mreg, 
+                       TRACKR.ptrack, 
+                       TRACKR.ctrack);
+           }
+           fluka->PrimaryIonisationStepping(nprim);
+       } // primary ionisation switched on
+    } // tracking resumed
 } // end of mgdraw
 } // end of extern "C"