Debug print removed.

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

diff --git a/TFluka/source.cxx b/TFluka/source.cxx
index 6c6f0142005a098d930133d71f408401e72406bd..0a8d907d7ca347db7e7ad94b9eb61bbb5804d966 100644 (file)
--- a/TFluka/source.cxx
+++ b/TFluka/source.cxx
@@ -5,14 +5,16 @@
 #include "Fdblprc.h"  //(DBLPRC) fluka common
 #include "Fdimpar.h"  //(DIMPAR) fluka parameters
 #include "Fsourcm.h"  //(EPISOR) fluka common
-#include "Fflkstk.h"  //(FLKSTK)  fluka common
-#include "Fsumcou.h"  //(SUMCOU)  fluka common
+#include "Fflkstk.h"  //(FLKSTK) fluka common
+#include "Fsumcou.h"  //(SUMCOU) fluka common
 #include "Fpaprop.h"  //(PAPROP) fluka common
 #include "Fltclcm.h"  //(LTCLCM) fluka common
 #include "Fopphst.h"  //(OPPHST) fluka common
-
+#include "Fioiocm.h"  //(IOIOCM) fluka common
+#include "Fbeamcm.h"  //(BEAMCM) fluka common
 //Virutal MC
 #include "TFluka.h"
+#include "TFlukaIon.h"
 
 #include "TVirtualMCStack.h"
 //#include "TVirtualMCApplication.h"
@@ -29,14 +31,19 @@
 # define georeg georeg_
 # define geohsm geohsm_
 # define soevsv soevsv_
+# define dcdion dcdion_
+# define setion setion_
+# define stisbm stisbm_
 #else
 # define source SOURCE
 # define geocrs GEOCRS
 # define georeg GEOREG
 # define geohsm GEOHSM
 # define soevsv SOEVSV
+# define dcdion DCDION
+# define setion SETION
+# define stisbm STISBM
 #endif
-
 extern "C" {
   //
   // Prototypes for FLUKA functions
@@ -46,6 +53,10 @@ extern "C" {
                            Int_t &, Int_t &);
   void type_of_call geohsm(Int_t &, Int_t &, Int_t &, Int_t &);
   void type_of_call soevsv();
+  void type_of_call dcdion(Int_t &);
+  void type_of_call setion(Int_t &);
+  void type_of_call stisbm(Int_t &, Int_t &, Int_t &);
+
  /*
    *----------------------------------------------------------------------*
    *                                                                      *
@@ -165,13 +176,58 @@ extern "C" {
     /* Cosines (tx,ty,tz)*/
     Double_t cosx = particle->Px()/particle->P();
     Double_t cosy = particle->Py()/particle->P();
-    Double_t cosz = TMath::Sqrt(oneone - cosx*cosx - cosy*cosy);
+    Double_t cosxy = cosx * cosx + cosy * cosy;
+    Double_t cosz;
+
+    if (cosxy < 1.) {
+       cosz = TMath::Sqrt(oneone - cosxy);
+    } else {
+       cosx /= TMath::Sqrt(cosxy);
+       cosy /= TMath::Sqrt(cosxy);     
+       cosz = 0.;
+    }
+    
+    
+    
     if (particle->Pz() < 0.) cosz = -cosz;    
 
     if (pdg != 50000050 &&  pdg !=  50000051) {
         FLKSTK.npflka++;
         Int_t ifl =  fluka-> IdFromPDG(pdg);
-        FLKSTK.iloflk[FLKSTK.npflka] = ifl;
+       Int_t ionid;
+       
+       if (ifl == -2) {
+           Int_t ia = TFlukaIon::GetA(pdg);
+           Int_t iz = TFlukaIon::GetZ(pdg);
+           Int_t is = TFlukaIon::GetIsomerNumber(pdg);
+           
+           if (is == 0) {
+               IOIOCM.iproa = ia;
+               IOIOCM.iproz = iz;
+               BEAMCM.ijhion = iz * 1000 + ia;
+               BEAMCM.ijhion = 100 * BEAMCM.ijhion + 30;
+               ionid = BEAMCM.ijhion;
+               dcdion(ionid);
+               setion(ionid);
+               FLKSTK.iloflk[FLKSTK.npflka] = ionid;
+               FLKSTK.lraddc[FLKSTK.npflka] = 0;
+           } else {
+               BEAMCM.lrdbea = 1;
+               stisbm(ia, iz, is);
+               BEAMCM.ijhion = iz * 1000 + ia;
+               BEAMCM.ijhion = 100 * BEAMCM.ijhion + 30;
+               ionid = BEAMCM.ijhion;
+               dcdion(ionid);
+               setion(ionid);
+               FLKSTK.iloflk[FLKSTK.npflka] = ionid;
+               FLKSTK.lraddc[FLKSTK.npflka] = 1;
+           }
+       } else {
+           FLKSTK.iloflk[FLKSTK.npflka] = ifl;
+           FLKSTK.lraddc[FLKSTK.npflka] = 0;
+           ionid = ifl;
+       }
+       
         /* Wtflk is the weight of the particle*/
         FLKSTK.wtflk[FLKSTK.npflka] = oneone;
         SUMCOU.weipri += FLKSTK.wtflk[FLKSTK.npflka];
@@ -207,7 +263,8 @@ extern "C" {
         
         /* Kinetic energy */
         Double_t p    = particle->P();
-        Double_t mass = PAPROP.am[ifl + 6];
+//        Double_t mass = PAPROP.am[ifl + 6];
+       Double_t mass = PAPROP.am[ionid + 6];
         FLKSTK.tkeflk[FLKSTK.npflka] =  TMath::Sqrt( p * p + mass * mass) - mass;
         /* Particle momentum*/
         FLKSTK.pmoflk [FLKSTK.npflka] = p;
@@ -323,6 +380,7 @@ extern "C" {
 //  Pre-track actions at for primary tracks
 //
     if (particleIsPrimary) {
+       fluka->SetCaller(kSODRAW);
         TVirtualMCApplication::Instance()->BeginPrimary();
         TVirtualMCApplication::Instance()->PreTrack();
     }