Support for heavy ions as primary in source.cxx
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 17 Apr 2008 21:56:07 +0000 (21:56 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 17 Apr 2008 21:56:07 +0000 (21:56 +0000)
TFluka/TFlukaIon.cxx
TFluka/TFlukaIon.h
TFluka/source.cxx

index c9cdb9f..11a4802 100644 (file)
@@ -62,6 +62,26 @@ Int_t TFlukaIon::GetIonPdg(Int_t z, Int_t a, Int_t i)
   return 1000000000 + 10*1000*z + 10*a + i;
 }  
 
+Int_t TFlukaIon::GetZ(Int_t pdg)
+{
+// Acording to
+// http://cepa.fnal.gov/psm/stdhep/pdg/montecarlorpp-2006.pdf
+
+  return (pdg - 1000000000)/10000; 
+}  
+
+
+Int_t TFlukaIon::GetA(Int_t pdg)
+{
+// Acording to
+// http://cepa.fnal.gov/psm/stdhep/pdg/montecarlorpp-2006.pdf
+
+    Int_t a = pdg - 1000000000;
+    a %= 10000;
+    a /= 10;
+    return (a);
+}  
+
 void TFlukaIon::AddIon(Int_t a, Int_t z)
 {
 
index ddc9ae2..7d1619f 100644 (file)
@@ -26,6 +26,7 @@ public:
     TFlukaIon(const char* name, Int_t z, Int_t a, Int_t q, Double_t exE, Double_t mass = 0.);
     Int_t    GetZ()                const  {return fZ;}
     Int_t    GetA()                const  {return fA;}
+
     Int_t    GetQ()                const  {return fQ;}
     Double_t GetExcitationEnergy() const  {return fExEnergy;}
     Double_t GetMass()             const  {return fMass;}
@@ -35,6 +36,9 @@ public:
     //
     static void  AddIon(Int_t a, Int_t z);
     static Int_t GetIonPdg(Int_t z, Int_t a, Int_t i = 0);
+    static Int_t    GetZ(Int_t pdg);
+    static Int_t    GetA(Int_t pdg);
+
  protected:
     Int_t    fZ;         // Z
     Int_t    fA;         // A
index 0e833e6..768e68c 100644 (file)
@@ -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"
 # define georeg georeg_
 # define geohsm geohsm_
 # define soevsv soevsv_
+# define dcdion dcdion_
+# define setion setion_
 #else
 # define source SOURCE
 # define geocrs GEOCRS
 # define georeg GEOREG
 # define geohsm GEOHSM
 # define soevsv SOEVSV
+# define dcdion DCDION
+# define setion SETION
 #endif
 
 extern "C" {
@@ -46,6 +52,9 @@ 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 &);
+
  /*
    *----------------------------------------------------------------------*
    *                                                                      *
@@ -183,7 +192,24 @@ extern "C" {
     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);
+           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;
+       } else {
+           FLKSTK.iloflk[FLKSTK.npflka] = ifl;
+           ionid = ifl;
+       }
+       
         /* Wtflk is the weight of the particle*/
         FLKSTK.wtflk[FLKSTK.npflka] = oneone;
         SUMCOU.weipri += FLKSTK.wtflk[FLKSTK.npflka];
@@ -219,7 +245,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;