- Improved support for heavy fragment transport
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 28 Nov 2007 15:46:20 +0000 (15:46 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 28 Nov 2007 15:46:20 +0000 (15:46 +0000)
- Improved IsAlive definition.

TFluka/TFluka.cxx
TFluka/TFluka.h
TFluka/mgdraw.cxx
TFluka/stuprf.cxx

index 20469f3..83756e7 100644 (file)
@@ -1849,7 +1849,10 @@ Double_t TFluka::Etot() const
               icode == kKASKADbrems || 
               icode == kKASKADpair) {
          return (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);      
+      } else {
+         return TRACKR.etrack;
       }
+      
   }
   else if (caller == kSODRAW) {
       Int_t ist  = FLKSTK.npflka;
@@ -1859,6 +1862,7 @@ Double_t TFluka::Etot() const
       Double_t e = TMath::Sqrt(p * p + m * m);
       return e;
   }
+  printf("Etot %5d %5d \n", caller, icode);
   
   return -1000.0;
 }
@@ -2119,7 +2123,6 @@ Int_t TFluka::StepProcesses(TArrayI &proc) const
   //
     FlukaProcessCode_t icode   = GetIcode();
     FlukaCallerCode_t  caller  = GetCaller();
-    
     proc.Set(1);
     TMCProcess iproc;
     if (caller == kBXEntering || caller == kBXExiting || caller == kEEDRAW || caller == kSODRAW) {
@@ -2165,6 +2168,16 @@ Int_t TFluka::StepProcesses(TArrayI &proc) const
        case kEMFSCOstopping2:
        case kKASNEUstopping:
            iproc = kPStop;
+           break; 
+       case kKASKADinelint:
+       case kKASNEUhadronic:
+           iproc = kPHadronic;
+           break;
+       case kKASKADinelarecoil:
+           iproc = kPHadronic;
+           break;
+       case kKASKADnelint:
+           iproc = kPHElastic;
            break;
        case kKASOPHabsorption:
            iproc = kPLightAbsorption;
@@ -2570,6 +2583,19 @@ void TFluka::AddParticlesToPdgDataBase() const
                      0,0,"Special",GetSpecialPdg(51));
 }
 
+void TFluka::AddIon(Int_t a, Int_t z) const
+{
+
+    // Add a new ion
+    TDatabasePDG *pdgDB = TDatabasePDG::Instance();
+    const Double_t kAu2Gev   = 0.9314943228;
+    Int_t pdg =  GetIonPdg(z, a);
+    if (pdgDB->GetParticle(pdg)) return;
+    
+    pdgDB->AddParticle(Form("Iion A  = %5d Z = %5d", a, z),"Ion", Float_t(a) * kAu2Gev + 8.071e-3, kTRUE,
+                      0, 3 * z, "Ion", pdg);
+}
+
 //
 // Info about primary ionization electrons
 //
index 5e4e545..746abdc 100644 (file)
@@ -400,6 +400,9 @@ class TFluka : public TVirtualMC {
   void     GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Double_t& z, Double_t& t) const;
   void     SetCurrentPrimaryElectronIndex(Int_t i)  {fPrimaryElectronIndex = i;}
   void     PrimaryIonisationStepping(Int_t nprim);
+
+  void  AddIon(Int_t a, Int_t z) const;
+  Int_t GetIonPdg(Int_t z, Int_t a, Int_t i = 0) const;
  private:
    
   // Copy constructor and operator= declared but not implemented (-Weff++ flag)
@@ -408,7 +411,7 @@ class TFluka : public TVirtualMC {
  
   void  PrintHeader();
   void  AddParticlesToPdgDataBase() const;
-  Int_t GetIonPdg(Int_t z, Int_t a, Int_t i = 0) const;
+
   Int_t GetSpecialPdg(Int_t number) const;
   
   Float_t* CreateFloatArray(Double_t* array, Int_t size) const;
index 6490139..1b57534 100644 (file)
@@ -45,7 +45,6 @@ void mgdraw(Int_t& icode, Int_t& mreg)
     }
     
     Int_t verbosityLevel = fluka->GetVerbosityLevel();
-
     if (TRACKR.jtrack < -6) {
        // from -7 to -12 = "heavy" fragment
        // assing parent id
@@ -54,7 +53,7 @@ void mgdraw(Int_t& icode, Int_t& mreg)
           cout << "mgdraw: (heavy fragment) jtrack < -6 =" << TRACKR.jtrack
                << " assign parent pdg=" << fluka->PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
        }
-       TRACKR.jtrack = TRACKR.ispusr[mkbmx2 - 3];
+//       TRACKR.jtrack = TRACKR.ispusr[mkbmx2 - 3];
     }
     
     cppstack->SetCurrentTrack(trackId);
index 3bd2e2a..a89eb92 100644 (file)
@@ -14,6 +14,7 @@
 #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
@@ -34,6 +35,8 @@ 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
@@ -62,8 +65,7 @@ extern "C" {
   FLKSTK.ispark[FLKSTK.npflka][mkbmx2 - 4] = TRACKR.ispusr[mkbmx2 - 1];  // current track number
   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;
   
@@ -86,14 +88,32 @@ extern "C" {
 
     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 pdg;
+    Int_t a, z;
+    
+    if (kpart < -6) {
+       kpart = -kpart;
+       z = kpart / 100000;
+       a = kpart - z * 100000;
+       a /= 100;
+       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 (kpart < 100000) {
+       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;
@@ -103,8 +123,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;
@@ -146,12 +164,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() 
             << "kin. energy [keV] = " <<  GENSTK.tki[numsec-1] * 1.e6 << endl << endl;
        
+    }
+    
 
 //
 //  Save current track number