Parent id and current track id handling corrected.

diff --git a/TFluka/stupre.cxx b/TFluka/stupre.cxx
index d4a582f60a8e22f11407b5ac8a6f84e90671ddfa..aec4c01c91649ce0d7b021b089f949c117a7b29d 100644 (file)
--- a/TFluka/stupre.cxx
+++ b/TFluka/stupre.cxx
@@ -67,13 +67,13 @@ void stupre()
 // Increment the track number and put it into the last flag
 
   Int_t kp;
-  for (kp=EMFSTK.npstrt-1; npnw<=EMFSTK.np-1; kp++) {
+  for (kp = EMFSTK.npstrt-1; npnw <= EMFSTK.np-1; kp++) {
 
 //* save the parent track number and reset it at each loop
-    Int_t ncrtrk = EVTFLG.ntrcks;
+    Int_t done = 1;
+
+    Int_t parent =  TRACKR.ispusr[mkbmx2-1];
     
-    Int_t done = 0;
-    Int_t parent = ncrtrk;
     Int_t flukaid = 0;
     if (EMFSTK.iq[kp] == -1) flukaid = 3;
     else if (EMFSTK.iq[kp] == 0)  flukaid = 7;
@@ -100,71 +100,71 @@ void stupre()
 //* all secondaries are true
     if ((EVTFLG.lpairp == 1) || (EVTFLG.lphoel == 1) ||
         (EVTFLG.lannfl == 1) || (EVTFLG.lannrs == 1)) {
-
-      if (EVTFLG.lpairp == 1) mech = kPPair;
-      else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
-      else mech = kPAnnihilation;
+       
+       if (EVTFLG.lpairp == 1) mech = kPPair;
+       else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
+       else mech = kPAnnihilation;
         cppstack->SetTrack(done, parent, pdg,
-           px, py, pz, e, vx, vy, vz, tof,
-           polx, poly, polz, mech, ntr, weight, is);
-
-cout << endl << " !!! stupre: ntr=" << ntr << endl;
-      EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+                          px, py, pz, e, vx, vy, vz, tof,
+                          polx, poly, polz, mech, ntr, weight, is);
+       
+       cout << endl << " !!! stupre: ntr=" << ntr  << " parent=" << parent << endl;
+       EMFSTK.iespak[kp][mkbmx2-1] = ntr;
     } // end of lpairp, lphoel, lannfl, lannrs
     
 //* Compton: secondary is true only if charged (e+, e-)
     else if ((EVTFLG.lcmptn == 1)) {
-      if (EMFSTK.iq[kp] != 0) {
-        mech = kPCompton;
-        cppstack->SetTrack(done, parent, pdg,
-           px, py, pz, e, vx, vy, vz, tof,
-           polx, poly, polz, mech, ntr, weight, is);
-cout << endl << " !!! stupre: ntr=" << ntr << endl;
-        EMFSTK.iespak[kp][mkbmx2-1] = ntr;
-      }
+       if (EMFSTK.iq[kp] != 0) {
+           mech = kPCompton;
+           cppstack->SetTrack(done, parent, pdg,
+                              px, py, pz, e, vx, vy, vz, tof,
+                              polx, poly, polz, mech, ntr, weight, is);
+           cout << endl << " !!! stupre: ntr=" << ntr  << " parent=" << parent << endl;
+           EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+       }
     } // end of lcmptn
-
+    
 //* Bremsstrahlung: true secondary only if charge = 0 (photon)
     else if ((EVTFLG.lbrmsp == 1)) {
-      if (EMFSTK.iq[kp] == 0) {
-        mech = kPBrem;
-        cppstack->SetTrack(done, parent, pdg,
-           px, py, pz, e, vx, vy, vz, tof,
-           polx, poly, polz, mech, ntr, weight, is);
-cout << endl << " !!! stupre: ntr=" << ntr << endl;
-        EMFSTK.iespak[kp][mkbmx2-1] = ntr;
-      }
+       if (EMFSTK.iq[kp] == 0) {
+           mech = kPBrem;
+           cppstack->SetTrack(done, parent, pdg,
+                              px, py, pz, e, vx, vy, vz, tof,
+                              polx, poly, polz, mech, ntr, weight, is);
+           cout << endl << " !!! stupre: ntr=" << ntr  << " parent=" << parent << endl;
+           EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+       }
     } // end of lbrmsp
-
+    
 //* Delta ray: If Bhabha, true secondary only if negative (electron)
     else if ((EVTFLG.ldltry == 1)) {
-      if (lbhabh == 1) {
-        if (EMFSTK.iq[kp] == 0) {
-          mech = kPDeltaRay;
-          cppstack->SetTrack(done, parent, pdg,
-            px, py, pz, e, vx, vy, vz, tof,
-            polx, poly, polz, mech, ntr, weight, is);
-          EMFSTK.iespak[kp][mkbmx2-1] = ntr;
-        } // end of Bhabha
-      }
-
+       if (lbhabh == 1) {
+           if (EMFSTK.iq[kp] == 0) {
+               mech = kPDeltaRay;
+               cppstack->SetTrack(done, parent, pdg,
+                                  px, py, pz, e, vx, vy, vz, tof,
+                                  polx, poly, polz, mech, ntr, weight, is);
+               EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+           } // end of Bhabha
+       }
+       
 //* Delta ray: Otherwise Moller: true secondary is the electron with
 //*            lower energy, which has been put higher in the stack
-      else if (kp == EMFSTK.np) {
-        mech = kPDeltaRay;
-        cppstack->SetTrack(done, parent, pdg,
-          px, py, pz, e, vx, vy, vz, tof,
-          polx, poly, polz, mech, ntr, weight, is);
-cout << endl << " !!! stupre: ntr=" << ntr << endl;
-        EMFSTK.iespak[kp][mkbmx2-1] = ntr;
-      } // end of Delta ray
+       else if (kp == EMFSTK.np) {
+           mech = kPDeltaRay;
+           cppstack->SetTrack(done, parent, pdg,
+                              px, py, pz, e, vx, vy, vz, tof,
+                              polx, poly, polz, mech, ntr, weight, is);
+           cout << endl << " !!! stupre: ntr=" << ntr << " parent=" << parent << endl;
+           EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+       } // end of Delta ray
     } // end of ldltry
-
+    
   } // end of loop
-
+  
 // !!! TO BE CONFIRMED !!!
   EVTFLG.ntrcks = EMFSTK.iespak[EMFSTK.np-1][mkbmx2-1];
-
+  
 } // end of stupre
 } // end of extern "C"
 

diff --git a/TFluka/stuprf.cxx b/TFluka/stuprf.cxx
index 49ab75b1c7d12bcd36260d7bf25e787afd7c715a..48898f28f19d5721709a2cdff95671a4e6f7ef61 100644 (file)
--- a/TFluka/stuprf.cxx
+++ b/TFluka/stuprf.cxx
@@ -1,5 +1,6 @@
 #include <Riostream.h>
 #include "AliRun.h"
+#include "AliStack.h"
 #include "TFluka.h"
 #ifndef WIN32
 # define stuprf stuprf_
@@ -47,10 +48,10 @@ void stuprf(Int_t& ij, Int_t& mreg,
 // TRACKR.spausr = user defined spare variables for the current particle
 // TRACKR.ispusr = user defined spare flags for the current particle
   Int_t ispr;
-  for (ispr=0; ispr<=mkbmx1-1; ispr++) {
+  for (ispr = 0; ispr <= mkbmx1 - 1; ispr++) {
     STACK.sparek[STACK.lstack][ispr] = TRACKR.spausr[ispr];
   }  
-  for (ispr=0; ispr<=mkbmx2-1; ispr++) {
+  for (ispr = 0; ispr <= mkbmx2 - 1; ispr++) {
     STACK.ispark[STACK.lstack][ispr] = TRACKR.ispusr[ispr];
   }  
  
@@ -61,52 +62,66 @@ void stuprf(Int_t& ij, Int_t& mreg,
   
 // EVTFLG.ntrcks = track number
 // Increment the track number and put it into the last flag
-  if (numsec-1 > npprmr) {
+// was numsec -1
+// clarify with Alberto
+  if (numsec > npprmr) {
 // Now call the SetTrack(...)
-    Int_t done = 0;
-    Int_t parent = TRACKR.ispusr[mkbmx2-1];
-    Int_t pdg = fluka->PDGFromId(ij);
-    
-    Double_t px = FINUC.plr[numsec-1]*FINUC.cxr[numsec-1];
-    Double_t pz = FINUC.plr[numsec-1]*FINUC.cyr[numsec-1];
-    Double_t py = FINUC.plr[numsec-1]*FINUC.czr[numsec-1];
+    Int_t done = 1;
+
+    Int_t parent =  TRACKR.ispusr[mkbmx2-1];
+
+    Int_t pdg = fluka->PDGFromId(FINUC.kpart[numsec-1]);
+    Double_t px = FINUC.plr[numsec-1] * FINUC.cxr[numsec-1];
+    Double_t pz = FINUC.plr[numsec-1] * FINUC.cyr[numsec-1];
+    Double_t py = FINUC.plr[numsec-1] * FINUC.czr[numsec-1];
     Double_t e  = FINUC.tki[numsec-1] + PAPROP.am[FINUC.kpart[numsec-1]+6];
     Double_t vx = xx;
     Double_t vy = yy;
     Double_t vz = zz;
-    Double_t tof = TRACKR.atrack;
+    
+    Double_t tof  = TRACKR.atrack;
     Double_t polx = FINUC.cxrpol[numsec-1];
     Double_t poly = FINUC.cyrpol[numsec-1];
     Double_t polz = FINUC.czrpol[numsec-1];
+    
 
     TMCProcess mech = kPHadronic;
-    if (EVTFLG.ldecay == 1) mech = kPDecay;
-    else if (EVTFLG.ldltry == 1) mech = kPDeltaRay;
-    else if (EVTFLG.lpairp == 1) mech = kPPair;
-    else if (EVTFLG.lbrmsp == 1) mech = kPBrem;
+    
+    if (EVTFLG.ldecay == 1) {
+       mech = kPDecay;
+       cout << endl << "Decay" << endl;
+       
+    } else if (EVTFLG.ldltry == 1) {
+       mech = kPDeltaRay;
+       cout << endl << "Delta Ray" << endl;
+       
+    } else if (EVTFLG.lpairp == 1) {
+       mech = kPPair;
+       cout << endl << "Pair Production" << endl;
+       
+    } else if (EVTFLG.lbrmsp == 1) {
+       mech = kPBrem;
+       cout << endl << "Bremsstrahlung" << endl;
+       
+    }
+    
+
     Double_t weight = FINUC.wei[numsec-1];
     Int_t is = 0;
     Int_t ntr;  
-
-//virtual void SetTrack(Int_t done, Int_t parent, Int_t pdg,
-//Double_t px, Double_t py, Double_t pz, Double_t e,
-//Double_t vx, Double_t vy, Double_t vz, Double_t tof,
-//Double_t polx, Double_t poly, Double_t polz,
-//TMCProcess mech, Int_t& ntr, Double_t weight,
-//Int_t is) = 0;
-
-    
+    // 
+    // Save particle in VMC stack
     cppstack->SetTrack(done, parent, pdg,
-                   px, py, pz, e,
-                   vx, vy, vz, tof,
-                   polx, poly, polz,
-                   mech, ntr, weight, is);
-
-cout << endl << " !!! stuprf: ntr=" << ntr << endl;
-    EVTFLG.ntrcks = ntr;
-    STACK.ispark[STACK.lstack][mkbmx2-1] = EVTFLG.ntrcks;
+                      px, py, pz, e,
+                      vx, vy, vz, tof,
+                      polx, poly, polz,
+                      mech, ntr, weight, is);
+    cout << endl << " !!! stuprf: ntr=" << ntr << "pdg " << pdg << " parent=" << parent << "numsec " 
+        << numsec << "npprmr " << npprmr << endl;
+//
+//  Save current track number
+    STACK.ispark[STACK.lstack][mkbmx2-1] = ntr;
   } // end of if (numsec-1 > npprmr)
-
 } // end of stuprf
 } // end of extern "C"