Correct loop over secondaries.

diff --git a/TFluka/stupre.cxx b/TFluka/stupre.cxx
index aec4c01c91649ce0d7b021b089f949c117a7b29d..259fb5f31115827f29d17085a0fde03f70ca104d 100644 (file)
--- a/TFluka/stupre.cxx
+++ b/TFluka/stupre.cxx
@@ -67,7 +67,7 @@ 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; kp <= EMFSTK.np - 1; kp++) {
 
 //* save the parent track number and reset it at each loop
     Int_t done = 1;
@@ -78,19 +78,19 @@ void stupre()
     if (EMFSTK.iq[kp] == -1) flukaid = 3;
     else if (EMFSTK.iq[kp] == 0)  flukaid = 7;
     else if (EMFSTK.iq[kp] == 0)  flukaid = 4;
-    Int_t pdg = fluka->PDGFromId(flukaid);
-    Double_t e  = EMFSTK.e[kp]*emvgev;
-    Double_t p = sqrt(e*e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
-    Double_t px = p * EMFSTK.u[kp];
-    Double_t pz = p * EMFSTK.v[kp];
-    Double_t py = p * EMFSTK.w[kp];
-    Double_t tof = EMFSTK.agemf[kp];
-    Double_t polx = EMFSTK.upol[kp];
-    Double_t poly = EMFSTK.vpol[kp];
-    Double_t polz = EMFSTK.wpol[kp];
-    Double_t vx = EMFSTK.x[kp];
-    Double_t vy = EMFSTK.y[kp];
-    Double_t vz = EMFSTK.z[kp];
+    Int_t pdg       = fluka->PDGFromId(flukaid);
+    Double_t e      = EMFSTK.e[kp] * emvgev;
+    Double_t p      = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
+    Double_t px     = p * EMFSTK.u[kp];
+    Double_t pz     = p * EMFSTK.v[kp];
+    Double_t py     = p * EMFSTK.w[kp];
+    Double_t tof    = EMFSTK.agemf[kp];
+    Double_t polx   = EMFSTK.upol[kp];
+    Double_t poly   = EMFSTK.vpol[kp];
+    Double_t polz   = EMFSTK.wpol[kp];
+    Double_t vx     = EMFSTK.x[kp];
+    Double_t vy     = EMFSTK.y[kp];
+    Double_t vz     = EMFSTK.z[kp];
     Double_t weight = EMFSTK.wt[kp];
     Int_t ntr;
     TMCProcess mech;
@@ -107,8 +107,8 @@ void stupre()
         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;
+       cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
+
        EMFSTK.iespak[kp][mkbmx2-1] = ntr;
     } // end of lpairp, lphoel, lannfl, lannrs
     
@@ -119,7 +119,7 @@ void stupre()
            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;
+           cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
            EMFSTK.iespak[kp][mkbmx2-1] = ntr;
        }
     } // end of lcmptn
@@ -131,7 +131,7 @@ void stupre()
            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;
+           cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
            EMFSTK.iespak[kp][mkbmx2-1] = ntr;
        }
     } // end of lbrmsp
@@ -139,23 +139,24 @@ void stupre()
 //* Delta ray: If Bhabha, true secondary only if negative (electron)
     else if ((EVTFLG.ldltry == 1)) {
        if (lbhabh == 1) {
-           if (EMFSTK.iq[kp] == 0) {
+           if (EMFSTK.iq[kp] == -1) {
                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;
+           cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
            } // end of Bhabha
-       }
+       } // lbhabh == 1
        
 //* 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) {
+       else if (kp == EMFSTK.np-1) {
            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;
+           cout << endl << " !!! stupre (Moller) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
            EMFSTK.iespak[kp][mkbmx2-1] = ntr;
        } // end of Delta ray
     } // end of ldltry
@@ -163,7 +164,7 @@ void stupre()
   } // end of loop
   
 // !!! TO BE CONFIRMED !!!
-  EVTFLG.ntrcks = EMFSTK.iespak[EMFSTK.np-1][mkbmx2-1];
+//  EVTFLG.ntrcks = EMFSTK.iespak[EMFSTK.np-1][mkbmx2-1];
   
 } // end of stupre
 } // end of extern "C"