Corrections for synchronisation of TFluka and FLUKA particle info during transport:

- gammas, e+/- below threshold
- case npprmr==1 in stuprf
- delta rays from ions
Remaining problems:
- ions with id , -6
- kaons changing cgarge
Ernesto Lopez Torres

diff --git a/TFluka/TFluka.cxx b/TFluka/TFluka.cxx
index 05661d221833545929e6fa295c1fbb9c317999ad..1f0d523b91c89be3194c027140f3a807b9203a7b 100644 (file)
--- a/TFluka/TFluka.cxx
+++ b/TFluka/TFluka.cxx
@@ -952,8 +952,8 @@ Int_t TFluka::PDGFromId(Int_t id) const
                printf("PDGFromId: Error intfluka < 0: %d\n", id);
            return -1;
        }
-       if (fVerbosityLevel >= 3)
-           printf("mpdgha called with %d %d \n", id, intfluka);
+//     if (fVerbosityLevel >= 3)
+//         printf("mpdgha called with %d %d \n", id, intfluka);
        return mpdgha(intfluka);
     } else {
        // ions and optical photons
@@ -1483,6 +1483,27 @@ Double_t TFluka::Edep() const
   }
 }
 
+//______________________________________________________________________________ 
+Int_t TFluka::CorrectFlukaId() const
+{
+   // since we don't put photons and e- created bellow transport cut on the vmc stack
+   // 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 also.
+   if( (IsTrackStop() )
+        && TRACKR.ispusr[mkbmx2 - 4] == TRACKR.ispusr[mkbmx2 - 1]
+        && TRACKR.jtrack != TRACKR.ispusr[mkbmx2 - 3] ) {
+      if (fVerbosityLevel >=3)
+         cout << "CorrectFlukaId() for icode=" << GetIcode()
+               << " track=" << TRACKR.ispusr[mkbmx2 - 1]
+               << " current PDG=" << PDGFromId(TRACKR.jtrack)
+               << " assign parent PDG=" << PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
+      return TRACKR.ispusr[mkbmx2 - 3]; // assign parent identity
+   }
+   return TRACKR.jtrack;
+}
+
+
 //______________________________________________________________________________ 
 Int_t TFluka::TrackPid() const
 {
@@ -1490,7 +1511,7 @@ Int_t TFluka::TrackPid() const
 // TRACKR.jtrack = identity number of the particle
   FlukaCallerCode_t caller = GetCaller();
   if (caller != kEEDRAW) {
-      return PDGFromId(TRACKR.jtrack);
+     return PDGFromId( CorrectFlukaId() );
   }
   else
     return -1000;
@@ -1504,7 +1525,7 @@ Double_t TFluka::TrackCharge() const
 // TRACKR.jtrack = identity number of the particle
   FlukaCallerCode_t caller = GetCaller();
   if (caller != kEEDRAW) 
-    return PAPROP.ichrge[TRACKR.jtrack+6];
+     return PAPROP.ichrge[CorrectFlukaId()+6];
   else
     return -1000.0;
 }
@@ -1516,7 +1537,7 @@ Double_t TFluka::TrackMass() const
 // TRACKR.jtrack = identity number of the particle
   FlukaCallerCode_t caller = GetCaller();
   if (caller != kEEDRAW)
-    return PAPROP.am[TRACKR.jtrack+6];
+     return PAPROP.am[CorrectFlukaId()+6];
   else
     return -1000.0;
 }
@@ -1634,15 +1655,15 @@ Bool_t   TFluka::IsTrackStop() const
 // True if the track energy has fallen below the threshold
 // means stopped by signal or below energy threshold
   FlukaProcessCode_t icode = GetIcode();
-  if (icode == kKASKADstopping  ||
-      icode == kKASKADtimekill  ||
-      icode == kEMFSCOstopping1 ||
-      icode == kEMFSCOstopping2 ||
-      icode == kEMFSCOtimekill  ||
-      icode == kKASNEUstopping  ||
-      icode == kKASNEUtimekill  ||
-      icode == kKASHEAtimekill  ||
-      icode == kKASOPHtimekill) return 1;
+  if (icode == kKASKADstopping  || // stopping particle
+      icode == kKASKADtimekill  || // time kill 
+      icode == kEMFSCOstopping1 || // below user-defined cut-off
+      icode == kEMFSCOstopping2 || // below user cut-off
+      icode == kEMFSCOtimekill  || // time kill
+      icode == kKASNEUstopping  || // neutron below threshold
+      icode == kKASNEUtimekill  || // time kill
+      icode == kKASHEAtimekill  || // time kill
+      icode == kKASOPHtimekill) return 1; // time kill
   else return 0;
 }
 

diff --git a/TFluka/TFluka.h b/TFluka/TFluka.h
index 8715af8e7f65f1e117ade3141f127f69791c59b9..de5c14d32d10b8f301038093d1b3b561816d1e51 100644 (file)
--- a/TFluka/TFluka.h
+++ b/TFluka/TFluka.h
@@ -222,6 +222,7 @@ class TFluka : public TVirtualMC {
   virtual Double_t TrackTime() const;
   virtual Double_t Edep() const;
   // Static properties
+  virtual Int_t    CorrectFlukaId() const;
   virtual Int_t    TrackPid() const;
   virtual Double_t TrackCharge() const;
   virtual Double_t TrackMass() const;

diff --git a/TFluka/bxdraw.cxx b/TFluka/bxdraw.cxx
index 56e211b0d1d31b5a8dee81b3748b03170d0d4afc..c9abc12bfd4604745265bfa0c79ff4be3f381fbe 100644 (file)
--- a/TFluka/bxdraw.cxx
+++ b/TFluka/bxdraw.cxx
@@ -34,7 +34,9 @@ void bxdraw(Int_t& icode, Int_t& mreg, Int_t& newreg,
        fluka->SetTrackIsExiting();
        fluka->SetCaller(kBXExiting);
        fluka->SetMreg(mreg,oldlttc);
-       (TVirtualMCApplication::Instance())->Stepping(); 
+       TVirtualMCStack* cppstack = fluka->GetStack();
+       cppstack->SetCurrentTrack( TRACKR.ispusr[mkbmx2-1] );
+       (TVirtualMCApplication::Instance())->Stepping();
     }
     if (newreg != fluka->GetDummyRegion()) {
        if (debug) printf("bxdraw (en) \n");
@@ -42,8 +44,10 @@ void bxdraw(Int_t& icode, Int_t& mreg, Int_t& newreg,
        fluka->SetTrackIsEntering();
        if (fluka->GetDummyBoundary() == 1) fluka->SetDummyBoundary(2);
        fluka->SetMreg(newreg,newlttc);
+       TVirtualMCStack* cppstack = fluka->GetStack();
+       cppstack->SetCurrentTrack( TRACKR.ispusr[mkbmx2-1] );
        (TVirtualMCApplication::Instance())->Stepping();
-    }      
+    }
 
 } // end of bxdraw
 } // end of extern "C"

diff --git a/TFluka/endraw.cxx b/TFluka/endraw.cxx
index 5a1610de37d246c707254239e439c15e45abdc81..a955075f35f248b6057e1f5e2efb2b74a064beaa 100644 (file)
--- a/TFluka/endraw.cxx
+++ b/TFluka/endraw.cxx
@@ -2,6 +2,7 @@
 #include "TVirtualMCApplication.h"
 #include "TGeoMaterial.h"
 #include "TGeoManager.h"
+#include <TParticle.h>
 #include "TFlukaCerenkov.h"
 
 #include "TFluka.h"
@@ -11,6 +12,7 @@
 #include "Ftrackr.h"  //(TRACKR) fluka common
 #include "Fltclcm.h"  //(LTCLCM) fluka common
 #include "Fpaprop.h"  //(PAPROP) fluka common
+
 #ifndef WIN32
 # define endraw endraw_
 #else
@@ -34,33 +36,49 @@ void endraw(Int_t& icode, Int_t& mreg, Double_t& rull, Double_t& xsco, Double_t&
   
   Float_t edep = rull;
   
-  if (icode == kKASKADinelarecoil) {
-    if (debug) cout << " For icode=" << icode << " Stepping is NOT called" << endl;
-    return;
-  }
-
   if (TRACKR.jtrack == -1) {
-// Handle quantum efficiency the G3 way
+  // Handle quantum efficiency the G3 way
       if (debug) printf("endraw: Cerenkov photon depositing energy: %d %e\n", mreg, rull);
       TGeoMaterial* material = (gGeoManager->GetCurrentVolume())->GetMaterial();
       TFlukaCerenkov*  cerenkov = dynamic_cast<TFlukaCerenkov*> (material->GetCerenkovProperties());
       if (cerenkov) {
-         Double_t eff = (cerenkov->GetQuantumEfficiency(rull));
-         if (gRandom->Rndm() > eff) {
-             edep = 0.;
-         }
+          Double_t eff = (cerenkov->GetQuantumEfficiency(rull));
+          if (gRandom->Rndm() > eff) {
+              edep = 0.;
+          }
       }
   }
 
+  TVirtualMCStack* cppstack = fluka->GetStack();
+  Int_t saveTrackId = cppstack->GetCurrentTrackNumber();
+
+  if (debug) {
+     cout << "ENDRAW For icode=" << icode << " stacktrack=" << saveTrackId
+          << " track=" << TRACKR.ispusr[mkbmx2-1] << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
+          << " edep="<< edep <<endl;
+  }
+
   if (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2) {
       fluka->SetIcode((FlukaProcessCode_t)icode);
       fluka->SetRull(edep);
+      if (icode == kKASKADelarecoil && TRACKR.ispusr[mkbmx2-5]) {
+         //  Elastic recoil and in stuprf npprmr > 0,
+         //  the secondary being loaded is actually still the interacting particle
+         cppstack->SetCurrentTrack( TRACKR.ispusr[mkbmx2-4] );
+         //      cout << "endraw elastic recoil track=" << TRACKR.ispusr[mkbmx2-1] << " parent=" << TRACKR.ispusr[mkbmx2-4]
+         //           << endl;
+      }
+      else
+         cppstack->SetCurrentTrack(TRACKR.ispusr[mkbmx2-1] );
       (TVirtualMCApplication::Instance())->Stepping();
+      
+//      cppstack->SetCurrentTrack( saveTrackId );
   } else {
   //
   // For icode 21,22 the particle has fallen below thresshold.
   // This has to be signalled to the StepManager() 
   //
+      cppstack->SetCurrentTrack( TRACKR.ispusr[mkbmx2-1] );
       fluka->SetRull(edep);
       fluka->SetIcode((FlukaProcessCode_t) icode);
       (TVirtualMCApplication::Instance())->Stepping();
@@ -68,6 +86,8 @@ void endraw(Int_t& icode, Int_t& mreg, Double_t& rull, Double_t& xsco, Double_t&
       fluka->SetIcode((FlukaProcessCode_t)icode);
       fluka->SetRull(0.);
       (TVirtualMCApplication::Instance())->Stepping();
+//      cppstack->SetCurrentTrack( saveTrackId );
+
   }
 } // end of endraw
 } // end of extern "C"

diff --git a/TFluka/mgdraw.cxx b/TFluka/mgdraw.cxx
index 221596b56640dfc8378e1a80662ba82b227473ff..5e89412aae5e29384699e26898531f56fe90ac58 100644 (file)
--- a/TFluka/mgdraw.cxx
+++ b/TFluka/mgdraw.cxx
@@ -24,7 +24,6 @@ void mgdraw(Int_t& icode, Int_t& mreg)
 {
     TFluka* fluka =  (TFluka*) gMC;
     if (mreg == fluka->GetDummyRegion()) return;
-    Int_t verbosityLevel = fluka->GetVerbosityLevel();
 //
 //  Make sure that stack has currrent track Id
 //
@@ -32,12 +31,25 @@ void mgdraw(Int_t& icode, Int_t& mreg)
     TVirtualMCStack* cppstack = fluka->GetStack();
     
     if (TRACKR.jtrack == -1) {
-       trackId = OPPHST.louopp[OPPHST.lstopp];
-       if (trackId == 0) {
-           trackId = FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1];
-       }
+        trackId = OPPHST.louopp[OPPHST.lstopp];
+        if (trackId == 0) {
+            trackId = FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1];
+        }
     } else {
-       trackId = TRACKR.ispusr[mkbmx2-1];
+        trackId = TRACKR.ispusr[mkbmx2-1];
+    }
+    
+    Int_t verbosityLevel = fluka->GetVerbosityLevel();
+
+    if (TRACKR.jtrack < -6) {
+       // (Unknow heavy Ion???)
+       // assing parent id ???
+       // id < -6 was skipped in stuprf =>   if (kpart < -6) return;
+       if (verbosityLevel >= 3) {
+          cout << "mgdraw: jtrack < -6 =" << TRACKR.jtrack
+               << " assign parent pdg=" << fluka->PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
+       }
+       TRACKR.jtrack = TRACKR.ispusr[mkbmx2 - 3];
     }
     
     cppstack->SetCurrentTrack(trackId);
@@ -50,38 +62,40 @@ void mgdraw(Int_t& icode, Int_t& mreg)
     fluka->SetCaller(kMGDRAW);
 
     if (!TRACKR.ispusr[mkbmx2 - 2]) {
-       //
-       // Single step
-       if (verbosityLevel >= 3) {
-           cout << endl << "mgdraw: energy deposition for:" << trackId << endl;
-       }
-       (TVirtualMCApplication::Instance())->Stepping();
-       fluka->SetTrackIsNew(kFALSE);
+        //
+        // Single step
+        if (verbosityLevel >= 3) {
+           cout << endl << "mgdraw: energy deposition for:" << trackId
+                 << " icode=" << icode
+                 << " pdg=" << fluka->PDGFromId(TRACKR.jtrack) << " flukaid="<< TRACKR.jtrack << endl;
+        }
+        (TVirtualMCApplication::Instance())->Stepping();
+        fluka->SetTrackIsNew(kFALSE);
     } else {
-       //
-       // Tracking is being resumed after secondary tracking
-       //
-       if (verbosityLevel >= 3) {
-           cout << endl << "mgdraw: resuming Stepping(): " << trackId << endl;
-       }
+        //
+        // Tracking is being resumed after secondary tracking
+        //
+        if (verbosityLevel >= 3) {
+            cout << endl << "mgdraw: resuming Stepping(): " << trackId << endl;
+        }
 
-       fluka->SetTrackIsNew(kTRUE);
-       fluka->SetCaller(kMGResumedTrack);
-       (TVirtualMCApplication::Instance())->Stepping();
+        fluka->SetTrackIsNew(kTRUE);
+        fluka->SetCaller(kMGResumedTrack);
+        (TVirtualMCApplication::Instance())->Stepping();
 
-       // Reset flag and stored values
-       TRACKR.ispusr[mkbmx2 - 2] = 0;
-       for (Int_t i = 0; i < 9; i++) TRACKR.spausr[i] = -1.;
+        // Reset flag and stored values
+        TRACKR.ispusr[mkbmx2 - 2] = 0;
+        for (Int_t i = 0; i < 9; i++) TRACKR.spausr[i] = -1.;
 
 
-       if (verbosityLevel >= 3) {
-           cout << endl << " !!! I am in mgdraw - first Stepping() after resume: " << icode << endl;
-           cout << endl << " Track Id = " << trackId << " region = " << mreg << endl;
-       }
+        if (verbosityLevel >= 3) {
+            cout << endl << " !!! I am in mgdraw - first Stepping() after resume: " << icode << endl;
+            cout << " Track= " << trackId << " region = " << mreg << endl;
+        }
 
-       fluka->SetTrackIsNew(kFALSE);
-       fluka->SetCaller(kMGDRAW);
-       (TVirtualMCApplication::Instance())->Stepping();
+        fluka->SetTrackIsNew(kFALSE);
+        fluka->SetCaller(kMGDRAW);
+        (TVirtualMCApplication::Instance())->Stepping();
     }
 } // end of mgdraw
 } // end of extern "C"

diff --git a/TFluka/stupre.cxx b/TFluka/stupre.cxx
index 937c2368ba718d1b5203bec9ac90b9a73de5d455..027497f5585350381c855533c1e2263117d79f31 100644 (file)
--- a/TFluka/stupre.cxx
+++ b/TFluka/stupre.cxx
@@ -92,6 +92,7 @@ void stupre()
     {
        EMFSTK.iespak[kp][mkbmx2-1] = TRACKR.ispusr[mkbmx2-1];
        EMFSTK.iespak[kp][mkbmx2-2] =  0;
+       EMFSTK.iespak[kp][mkbmx2-3] = TRACKR.jtrack;
        continue;
     }
 

diff --git a/TFluka/stuprf.cxx b/TFluka/stuprf.cxx
index 8eb2cdd4fe963680267429d16173986e9c0387ed..e80f2d1c6b32e26f8879910da570366316140606 100644 (file)
--- a/TFluka/stuprf.cxx
+++ b/TFluka/stuprf.cxx
@@ -15,6 +15,7 @@
 #include "Ftrackr.h"  //(TRACKR) fluka common
 #include "Fgenstk.h"  //(GENSTK)  fluka common
 
+
 //Virtual MC
 #include "TFluka.h"
 
@@ -25,8 +26,8 @@
 
 extern "C" {
     void stuprf(Int_t& /*ij*/, Int_t& /*mreg*/,
-               Double_t& xx, Double_t& yy, Double_t& zz,
-               Int_t& numsec, Int_t& npprmr)
+                Double_t& xx, Double_t& yy, Double_t& zz,
+                Int_t& numsec, Int_t& npprmr)
 {
 //*----------------------------------------------------------------------*
 //*                                                                      *
@@ -53,6 +54,11 @@ extern "C" {
     FLKSTK.ispark[FLKSTK.npflka][ispr] = TRACKR.ispusr[ispr];
   }  
  
+  // save parent info
+  FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 3] = TRACKR.jtrack;   // fluka particle id
+  FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 4] = TRACKR.ispusr[mkbmx2 - 1];  // current track number
+  FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 5] = npprmr; // flag special case when npprmr>0
+
 // Get the pointer to the VMC
   TFluka* fluka =  (TFluka*) gMC;
   Int_t verbosityLevel = fluka->GetVerbosityLevel();
@@ -66,7 +72,11 @@ extern "C" {
 // Increment the track number and put it into the last flag
 // was numsec -1
 // clarify with Alberto
-  if (numsec > npprmr) {
+  
+//  npprmr > 0, the secondary being loaded is actually still the interacting
+//  particle (it can happen in some biasing situations)
+
+  if (numsec > npprmr || npprmr > 0) {
 // Now call the PushTrack(...)
     Int_t done = 0;
 
@@ -76,7 +86,6 @@ extern "C" {
 
     Int_t 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];
@@ -85,51 +94,71 @@ extern "C" {
     Double_t vx = xx;
     Double_t vy = yy;
     Double_t vz = zz;
-    
+
     Double_t tof  = TRACKR.atrack;
     Double_t polx = GENSTK.cxrpol[numsec-1];
     Double_t poly = GENSTK.cyrpol[numsec-1];
     Double_t polz = GENSTK.czrpol[numsec-1];
-    
+
 
     TMCProcess mech = kPHadronic;
-    
+
     if (EVTFLG.ldecay == 1) {
-       mech = kPDecay;
-       if (debug) cout << endl << "Decay" << endl;
-       
+        mech = kPDecay;
+        if (debug) cout << endl << "Decay" << endl;
     } else if (EVTFLG.ldltry == 1) {
-       mech = kPDeltaRay;
-       if (debug) cout << endl << "Delta Ray" << endl;
-       
+        mech = kPDeltaRay;
+        if( fluka->GetIcode() == kKASHEA ) {
+           //  For all interactions secondaries are put on GENSTK common (kp=1,np)
+           //  but for KASHEA delta ray generation where only the secondary elec-
+           //  tron is present and stacked on FLKSTK common for kp=lstack
+           pdg  = fluka->PDGFromId( FLKSTK.iloflk[FLKSTK.npflka] );
+           px   = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.txflk[FLKSTK.npflka];
+           py   = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.tyflk[FLKSTK.npflka];
+           pz   = FLKSTK.pmoflk[FLKSTK.npflka] * FLKSTK.tzflk[FLKSTK.npflka];
+           e    = FLKSTK.tkeflk[FLKSTK.npflka] + PAPROP.am[FLKSTK.iloflk[FLKSTK.npflka]+6];
+           polx = FLKSTK.txpol[FLKSTK.npflka];
+           poly = FLKSTK.typol[FLKSTK.npflka];
+           polz = FLKSTK.tzpol[FLKSTK.npflka];
+           if (debug) cout << endl << "Delta Ray from KASHEA...." << " pdg from FLKSTK=" << pdg << endl;
+        } else
+           if (debug) cout << endl << "Delta Ray" << endl;
     } else if (EVTFLG.lpairp == 1) {
-       mech = kPPair;
-       if (debug) cout << endl << "Pair Production" << endl;
-       
+        mech = kPPair;
+        if (debug) cout << endl << "Pair Production" << endl;
     } else if (EVTFLG.lbrmsp == 1) {
-       mech = kPBrem;
-       if (debug) cout << endl << "Bremsstrahlung" << endl;
-       
+        mech = kPBrem;
+        if (debug) cout << endl << "Bremsstrahlung" << endl;
     }
-    
 
     Double_t weight = GENSTK.wei[numsec-1];
     Int_t is = 0;
-    Int_t ntr;  
+    Int_t ntr;
     // 
     // Save particle in VMC stack
     cppstack->PushTrack(done, parent, pdg,
-                      px, py, pz, e,
-                      vx, vy, vz, tof,
-                      polx, poly, polz,
-                      mech, ntr, weight, is);
-    if (debug) cout << endl << " !!! stuprf: ntr=" << ntr << "pdg " << pdg << " parent=" << parent << "numsec " 
-        << numsec << "npprmr " << npprmr << endl;
+                       px, py, pz, e,
+                       vx, vy, vz, tof,
+                       polx, poly, polz,
+                       mech, ntr, weight, is);
+    if (debug)
+       cout << endl << " !!! stuprf: ntr=" << ntr << " pdg " << pdg << " parent=" << parent
+             << " parent_pdg="<< fluka->PDGFromId(TRACKR.jtrack) << " numsec "
+             << numsec << " npprmr " << npprmr << " icode=" << fluka->GetIcode()
+             << endl;
+
 //
 //  Save current track number
     FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1] = ntr;
     FLKSTK.ispark[FLKSTK.npflka][mkbmx2-2] = 0;
-  } // end of if (numsec-1 > npprmr)
+  } // end of if (numsec > npprmr)
+//  else {
+//     if(debug) {
+//        cout << endl << " !!! stuprf: skipping pushtrack   track=" << TRACKR.ispusr[mkbmx2-1]
+//              << " pdg " << fluka->PDGFromId(TRACKR.jtrack) << " numsec=" << numsec<< " npprmr=" << npprmr
+//              << " GENSTK pdg=" << fluka->PDGFromId(GENSTK.kpart[numsec-1]) << endl;
+//     }
+//  }
 } // end of stuprf
 } // end of extern "C"
 

diff --git a/TFluka/usdraw.cxx b/TFluka/usdraw.cxx
index ea9422db75e527818d2e712d0ef7a1e631c9d389..3074d480d427fe12fee7e4bdfbaf3c0a6519d932 100644 (file)
--- a/TFluka/usdraw.cxx
+++ b/TFluka/usdraw.cxx
@@ -43,8 +43,6 @@ void usdraw(Int_t& icode, Int_t& mreg,
          } // Track found in stack
       } // Loop over emf stack 
   } // Electromagnetic process
-  
-
 
   Int_t mlttc = LTCLCM.mlatm1;
   fluka->SetMreg(mreg, mlttc);
@@ -52,11 +50,16 @@ void usdraw(Int_t& icode, Int_t& mreg,
   fluka->SetYsco(ysco);
   fluka->SetZsco(zsco);
 
-  if (debug) printf("USDRAW: Number of track segments:%6d %6d %6d %10.3e\n", TRACKR.ntrack, TRACKR.mtrack, icode, TRACKR.atrack);
+  if (debug) printf("USDRAW: Number of track segments:%6d %6d icode=%d tof=%10.3e track=%d pdg=%d\n",
+  TRACKR.ntrack, TRACKR.mtrack, icode, TRACKR.atrack, TRACKR.ispusr[mkbmx2-1], fluka->PDGFromId(TRACKR.jtrack) );
+
+  TVirtualMCStack* cppstack = fluka->GetStack();
+  cppstack->SetCurrentTrack( TRACKR.ispusr[mkbmx2-1] );
 
   (TVirtualMCApplication::Instance())->Stepping();
   fluka->SetTrackIsNew(kFALSE);
-  
+
+ 
 } // end of usdraw
 } // end of extern "C"