Updated misalignment macros (Raffaele)

[u/mrichter/AliRoot.git] / TFluka / stuprf.cxx

diff --git a/TFluka/stuprf.cxx b/TFluka/stuprf.cxx
index 3e73b7a22d742ad300135791b8e2f52e516b3626..ae11875bead3c8d0c1d294c37619560beb6e587c 100644 (file)
--- a/TFluka/stuprf.cxx
+++ b/TFluka/stuprf.cxx
@@ -1,8 +1,11 @@
 #include <Riostream.h>
+#include "TCallf77.h"      //For the fortran calls
 #ifndef WIN32
 # define stuprf stuprf_
+# define usrdci usrdci_
 #else
 # define stuprf STUPRF
+# define usrdci USRDCI
 #endif
 //
 // Fluka include
@@ -14,17 +17,20 @@
 #include "Fflkstk.h"  //(FLKSTK)  fluka common
 #include "Ftrackr.h"  //(TRACKR) fluka common
 #include "Fgenstk.h"  //(GENSTK)  fluka common
+#include "Ffheavy.h"  //(FHEAVY)  fluka common
 
 
 //Virtual MC
 #include "TFluka.h"
-
 #include "TVirtualMCStack.h"
 #include "TVirtualMCApplication.h"
 #include "TParticle.h"
 #include "TVector3.h"
 
 extern "C" {
+
+    void type_of_call usrdci(int&, int&, int&, int&);
+
     void stuprf(Int_t& /*ij*/, Int_t& /*mreg*/,
                 Double_t& xx, Double_t& yy, Double_t& zz,
                 Int_t& numsec, Int_t& npprmr)
@@ -34,11 +40,14 @@ extern "C" {
 //*  SeT User PRoperties for Fluka particles                             *
 //*                                                                      *
 //*----------------------------------------------------------------------*
+// Get the pointer to the VMC
+  TFluka* fluka =  (TFluka*) gMC;
 
 // FLKSTK.npflka  = stack pointer
 // FLKSTK.louse   = user flag
 // TRACKR.llouse = user defined flag for the current particle
-  FLKSTK.louse[FLKSTK.npflka] = TRACKR.llouse;
+    
+   FLKSTK.louse[FLKSTK.npflka] = TRACKR.llouse;
 
 // mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR
 // mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR
@@ -47,22 +56,23 @@ extern "C" {
 // 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++) {
-    FLKSTK.sparek[FLKSTK.npflka][ispr] = TRACKR.spausr[ispr];
-  }  
-  for (ispr = 0; ispr <= mkbmx2 - 1; ispr++) {
-    FLKSTK.ispark[FLKSTK.npflka][ispr] = TRACKR.ispusr[ispr];
-  }  
- 
+  if (numsec <= npprmr) {
+      for (ispr = 0; ispr <= mkbmx1 - 1; ispr++) {
+         FLKSTK.sparek[FLKSTK.npflka][ispr] = TRACKR.spausr[ispr];
+      }  
+      for (ispr = 0; ispr <= mkbmx2 - 1; ispr++) {
+         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 - 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
+  FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 5] = npprmr;                     // flag npprmr>0
+
 
-// Get the pointer to the VMC
-  TFluka* fluka =  (TFluka*) gMC;
   Int_t verbosityLevel = fluka->GetVerbosityLevel();
-  Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
+  Bool_t debug = (verbosityLevel>=3)? kTRUE:kFALSE;
   
   fluka->SetTrackIsNew(kTRUE);
 //  TVirtualMC* fluka = TFluka::GetMC();
@@ -81,16 +91,36 @@ extern "C" {
 // Now call the PushTrack(...)
     Int_t done = 0;
 
-    Int_t parent =  TRACKR.ispusr[mkbmx2-1];
-    Int_t kpart  = GENSTK.kpart[numsec-1];
-    if (kpart < -6) return; // -7 to -12 = "heavy" fragment
-    Int_t  pdg  = fluka->PDGFromId(kpart);
-     
+    Int_t parent  = TRACKR.ispusr[mkbmx2-1];
+    Int_t kpart   = GENSTK.kpart[numsec-1];
+    Int_t pdg;
+    Int_t a = 1;
+    Int_t z = 1;
+    Int_t ism = 0;
+    if (kpart < -6) {
+//     kpart = -kpart;
+//     z = kpart / 100000;
+//     a = kpart - z * 100000;
+//     a /= 100;
+       usrdci(kpart, a, z, ism);
+       pdg = fluka->GetIonPdg(z, a);
+       fluka->AddIon(a, z);
+    } else {
+       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];
-    Double_t e  = GENSTK.tki[numsec-1] + PAPROP.am[GENSTK.kpart[numsec-1]+6];
-
+    Double_t e;
+    
+    if (a == 1) {
+       e  = GENSTK.tki[numsec-1] + PAPROP.am[GENSTK.kpart[numsec-1]+6];
+    } else {
+       e  = GENSTK.tki[numsec-1] + Double_t(a) * 0.93149 + 8.071e-3;
+    }
+    
+    
     Double_t vx = xx;
     Double_t vy = yy;
     Double_t vz = zz;
@@ -100,7 +130,6 @@ extern "C" {
     Double_t poly = GENSTK.cyrpol[numsec-1];
     Double_t polz = GENSTK.czrpol[numsec-1];
     
-
     TMCProcess mech = kPHadronic;
     if (EVTFLG.ldecay == 1) {
         mech = kPDecay;
@@ -142,12 +171,15 @@ extern "C" {
                        vx, vy, vz, tof,
                        polx, poly, polz,
                        mech, ntr, weight, is);
-    if (debug)
+    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
-            << endl;
+            << numsec << " npprmr " << npprmr << " icode=" << fluka->GetIcode() 
+            << "kin. energy [keV] = " <<  GENSTK.tki[numsec-1] * 1.e6 << endl << endl;
        
+    }
+    
 
 //
 //  Save current track number
@@ -162,5 +194,6 @@ extern "C" {
 //     }
 //  }
 } // end of stuprf
+
 } // end of extern "C"