]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TFluka/TFluka.cxx
Cerenkov photon update.
[u/mrichter/AliRoot.git] / TFluka / TFluka.cxx
index 1f27bcb2bef73c342d7ec6e030acc5994f1731cf..f9808f00a2acdae803ce2e5daa98c2e7929aff6e 100644 (file)
@@ -133,7 +133,6 @@ TFluka::TFluka()
    fTrackIsExiting(kFALSE),
    fTrackIsNew(kFALSE),
    fFieldFlag(kTRUE),
-   fGeneratePemf(kFALSE),
    fDummyBoundary(kFALSE),
    fStopped(kFALSE),
    fStopEvent(kFALSE),
@@ -151,6 +150,7 @@ TFluka::TFluka()
   //
   // Default constructor
   //
+    for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
 } 
  
 //______________________________________________________________________________ 
@@ -171,7 +171,6 @@ TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupporte
    fTrackIsExiting(kFALSE),
    fTrackIsNew(kFALSE),
    fFieldFlag(kTRUE),
-   fGeneratePemf(kFALSE),
    fDummyBoundary(kFALSE),
    fStopped(kFALSE),
    fStopEvent(kFALSE),
@@ -187,6 +186,8 @@ TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupporte
    fUserScore(new TObjArray(100)) 
 {
   // create geometry interface
+    for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
+    
    if (fVerbosityLevel >=3)
        cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
    SetCoreInputFileName();
@@ -205,8 +206,8 @@ TFluka::~TFluka()
         cout << "<== TFluka::~TFluka() destructor called." << endl;
     if (fMaterials) delete [] fMaterials;
     
-    delete fGeom;
-    delete fMCGeo;
+//    delete fGeom;
+//    delete fMCGeo;
     
     if (fUserConfig) {
         fUserConfig->Delete();
@@ -297,13 +298,7 @@ void TFluka::BuildPhysics() {
         }
     }
     
-    //
-    // At this stage we have the information on materials and cuts available.
-    // Now create the pemf file
-    
-    if (fGeneratePemf) fGeom->CreatePemfFile();
-    
-    //
+
     // Prepare input file with the current physics settings
     
     InitPhysics(); 
@@ -354,7 +349,7 @@ Bool_t TFluka::ProcessRun(Int_t nevent) {
 //
 // Run steering
 //
-
+    
   if (fVerbosityLevel >=3)
     cout << "==> TFluka::ProcessRun(" << nevent << ") called." 
          << endl;
@@ -500,6 +495,8 @@ void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
   // In this case, WMAT in output is changed to relative
   // weigths.
   //
+    printf("Mixture %5d %10s %5d \n", kmat, name, nlmat);
+    
   Int_t i,j;
   if (nlmat < 0) {
      nlmat = - nlmat;
@@ -596,7 +593,8 @@ void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
      delete [] wmatnew;
      return;     
    }
-   gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
+  printf("Mixture (2) %5d %10s %5d \n", kmat, name, nlmat);
+  gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
 } 
 
 //______________________________________________________________________________ 
@@ -885,7 +883,8 @@ void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
 //
 //  
 //  Create object holding Cerenkov properties
-//  
+// 
+    
     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
 //
 //  Pass object to medium
@@ -917,16 +916,61 @@ void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
 
 
 //______________________________________________________________________________ 
-void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/,
-                         Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
+void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
+                         Double_t *absco, Double_t *effic, Double_t *rindex) {
+//
+// Set Cerenkov properties for medium itmed
+//
+// npckov: number of sampling points
+// ppckov: energy values
+// absco:  absorption length
+// effic:  quantum efficiency
+// rindex: refraction index
 //
-//  Double_t version not implemented
+
+//
+//  Double_t version 
+  Float_t* fppckov = CreateFloatArray(ppckov, npckov);
+  Float_t* fabsco  = CreateFloatArray(absco,  npckov);
+  Float_t* feffic  = CreateFloatArray(effic,  npckov);
+  Float_t* frindex = CreateFloatArray(rindex, npckov);
+
+  SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex);
+
+  delete [] fppckov;
+  delete [] fabsco;
+  delete [] feffic;
+  delete [] frindex;
 }  
 
-void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t* /*ppckov*/,
-                         Double_t* /*absco*/, Double_t* /*effic*/, Double_t* /*rindex*/, Double_t* /*rfl*/) {
+void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t* ppckov,
+                         Double_t* absco, Double_t* effic, Double_t* rindex, Double_t* rfl) {
+//
+// Set Cerenkov properties for medium itmed
+//
+// npckov: number of sampling points
+// ppckov: energy values
+// absco:  absorption length
+// effic:  quantum efficiency
+// rindex: refraction index
+// rfl:    reflectivity for boundary to medium itmed
+//
+
 //
-// //  Double_t version not implemented
+// //  Double_t version 
+  Float_t* fppckov = CreateFloatArray(ppckov, npckov);
+  Float_t* fabsco  = CreateFloatArray(absco,  npckov);
+  Float_t* feffic  = CreateFloatArray(effic,  npckov);
+  Float_t* frindex = CreateFloatArray(rindex, npckov);
+  Float_t* frfl    = CreateFloatArray(rfl,    npckov);
+
+  SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex, frfl);
+
+  delete [] fppckov;
+  delete [] fabsco;
+  delete [] feffic;
+  delete [] frindex;
+  delete [] frfl;
 }
 
 // Euclid
@@ -948,7 +992,11 @@ Int_t TFluka::GetMedium() const {
 //
 //  Get the medium number for the current fluka region
 //
-    return fGeom->GetMedium(); // this I need to check due to remapping !!!
+    if (gGeoManager->IsOutside()) {
+       return (-1);
+    } else {
+       return (fGeom->GetMedium()); // this I need to check due to remapping !!!
+    }
 }
 
 //____________________________________________________________________________ 
@@ -1336,10 +1384,10 @@ void TFluka::TrackPosition(TLorentzVector& position) const
   if (caller == kENDRAW    || caller == kUSDRAW || 
       caller == kBXExiting || caller == kBXEntering || 
       caller == kUSTCKV) { 
-    position.SetX(GetXsco());
-    position.SetY(GetYsco());
-    position.SetZ(GetZsco());
-    position.SetT(TRACKR.atrack);
+      position.SetX(GetXsco());
+      position.SetY(GetYsco());
+      position.SetZ(GetZsco());
+      position.SetT(TRACKR.atrack);
   }
   else if (caller == kMGDRAW) {
       Int_t i = -1;
@@ -1359,18 +1407,19 @@ void TFluka::TrackPosition(TLorentzVector& position) const
       }
   }
   else if (caller == kSODRAW) { 
-    position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
-    position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
-    position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
-    position.SetT(0);
+      Int_t ist = FLKSTK.npflka;
+      position.SetX(FLKSTK.xflk[ist]);
+      position.SetY(FLKSTK.yflk[ist]);
+      position.SetZ(FLKSTK.zflk[ist]);
+      position.SetT(FLKSTK.agestk[ist]);
   } else if (caller == kMGResumedTrack) { 
-    position.SetX(TRACKR.spausr[0]);
-    position.SetY(TRACKR.spausr[1]);
-    position.SetZ(TRACKR.spausr[2]);
-    position.SetT(TRACKR.spausr[3]);
+      position.SetX(TRACKR.spausr[0]);
+      position.SetY(TRACKR.spausr[1]);
+      position.SetZ(TRACKR.spausr[2]);
+      position.SetT(TRACKR.spausr[3]);
   }
   else
-    Warning("TrackPosition","position not available");
+      Warning("TrackPosition","position not available");
 }
 
 //______________________________________________________________________________ 
@@ -1390,7 +1439,7 @@ void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
       y = GetYsco();
       z = GetZsco();
   }
-  else if (caller == kMGDRAW || caller == kSODRAW) { 
+  else if (caller == kMGDRAW) { 
       Int_t i = -1;
       if ((i = fPrimaryElectronIndex) > -1) {
          GetPrimaryElectronPosition(i, x, y, z);
@@ -1400,13 +1449,19 @@ void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
          z = TRACKR.ztrack[TRACKR.ntrack];
       }
   }
+  else if (caller == kSODRAW) { 
+      Int_t ist = FLKSTK.npflka;
+      x = FLKSTK.xflk[ist];
+      y = FLKSTK.yflk[ist];
+      z = FLKSTK.zflk[ist];
+  }
   else if (caller == kMGResumedTrack) {
-    x = TRACKR.spausr[0];
-    y = TRACKR.spausr[1];
-    z = TRACKR.spausr[2];
+      x = TRACKR.spausr[0];
+      y = TRACKR.spausr[1];
+      z = TRACKR.spausr[2];
   }
   else
-    Warning("TrackPosition","position not available");
+      Warning("TrackPosition","position not available");
 }
 
 //______________________________________________________________________________ 
@@ -1423,34 +1478,74 @@ void TFluka::TrackMomentum(TLorentzVector& momentum) const
   FlukaCallerCode_t  caller = GetCaller();
   FlukaProcessCode_t icode  = GetIcode();
   
-  if (caller != kEEDRAW && caller != kMGResumedTrack && 
+  if (caller  != kEEDRAW         && 
+      caller  != kMGResumedTrack && 
+      caller  != kSODRAW         &&
+      caller  != kUSDRAW         &&
       (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
-    if (TRACKR.ptrack >= 0) {
-      momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
-      momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
-      momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
-      momentum.SetE(TRACKR.etrack);
-      return;
-    }
-    else {
-      Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
-      momentum.SetPx(p*TRACKR.cxtrck);
-      momentum.SetPy(p*TRACKR.cytrck);
-      momentum.SetPz(p*TRACKR.cztrck);
-      momentum.SetE(TRACKR.etrack);
-      return;
-    }
+      if (TRACKR.ptrack >= 0) {
+         momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
+         momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
+         momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
+         momentum.SetE(TRACKR.etrack);
+         return;
+      }
+      else {
+         Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
+         momentum.SetPx(p*TRACKR.cxtrck);
+         momentum.SetPy(p*TRACKR.cytrck);
+         momentum.SetPz(p*TRACKR.cztrck);
+         momentum.SetE(TRACKR.etrack);
+         return;
+      }
   } else if  (caller == kMGResumedTrack) {
-    momentum.SetPx(TRACKR.spausr[4]);
-    momentum.SetPy(TRACKR.spausr[5]);
-    momentum.SetPz(TRACKR.spausr[6]);
-    momentum.SetE (TRACKR.spausr[7]);
-    return;
+      momentum.SetPx(TRACKR.spausr[4]);
+      momentum.SetPy(TRACKR.spausr[5]);
+      momentum.SetPz(TRACKR.spausr[6]);
+      momentum.SetE (TRACKR.spausr[7]);
+      return;
   } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
       momentum.SetPx(0.);
       momentum.SetPy(0.);
       momentum.SetPz(0.);
       momentum.SetE(TrackMass());
+      
+  } else if (caller == kSODRAW) {
+      Int_t ist  = FLKSTK.npflka;
+      Double_t p = FLKSTK.pmoflk[ist];
+      Int_t ifl  = FLKSTK.iloflk[ist];
+      Double_t m = PAPROP.am[ifl + 6];
+      Double_t e = TMath::Sqrt(p * p + m * m);
+      momentum.SetPx(p * FLKSTK.txflk[ist]);
+      momentum.SetPy(p * FLKSTK.tyflk[ist]);
+      momentum.SetPz(p * FLKSTK.tzflk[ist]);
+      momentum.SetE(e);
+  } else if (caller == kUSDRAW) {
+      if (icode == kEMFSCObrems  || 
+         icode == kEMFSCOmoller || 
+         icode == kEMFSCObhabha || 
+         icode == kEMFSCOcompton ) 
+      {
+         momentum.SetPx(fPint[0]);
+         momentum.SetPy(fPint[1]);
+         momentum.SetPz(fPint[2]);
+         momentum.SetE(fPint[3]);
+      } else if (icode == kKASKADdray  || 
+                icode == kKASKADbrems || 
+                icode == kKASKADpair) {
+         momentum.SetPx(GENSTK.plr[0] * GENSTK.cxr[0]);
+         momentum.SetPy(GENSTK.plr[0] * GENSTK.cyr[0]);
+         momentum.SetPz(GENSTK.plr[0] * GENSTK.czr[0]);
+         momentum.SetE (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);
+      } else {
+         Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack 
+                           - ParticleMassFPC(TRACKR.jtrack) 
+                           * ParticleMassFPC(TRACKR.jtrack));
+         momentum.SetPx(p*TRACKR.cxtrck);
+         momentum.SetPy(p*TRACKR.cytrck);
+         momentum.SetPz(p*TRACKR.cztrck);
+         momentum.SetE(TRACKR.etrack);
+      }
   }
   else
     Warning("TrackMomentum","momentum not available");
@@ -1469,7 +1564,10 @@ void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e
 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
   FlukaCallerCode_t   caller = GetCaller();
   FlukaProcessCode_t  icode  = GetIcode();
-  if (caller != kEEDRAW && caller != kMGResumedTrack && 
+  if (caller != kEEDRAW         && 
+      caller != kMGResumedTrack && 
+      caller != kSODRAW         &&
+      caller != kUSDRAW         &&
       (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
     if (TRACKR.ptrack >= 0) {
       px = TRACKR.ptrack*TRACKR.cxtrck;
@@ -1497,9 +1595,42 @@ void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e
       py = 0.;
       pz = 0.;
       e  = TrackMass();
+  } else if (caller == kSODRAW) {
+      Int_t ist  = FLKSTK.npflka;
+      Double_t p = FLKSTK.pmoflk[ist];
+      Int_t ifl  = FLKSTK.iloflk[ist];
+      Double_t m = PAPROP.am[ifl + 6];
+               e = TMath::Sqrt(p * p + m * m);
+      px = p * FLKSTK.txflk[ist];
+      py = p * FLKSTK.tyflk[ist];
+      pz = p * FLKSTK.tzflk[ist];
+  } else if (caller == kUSDRAW) {
+      if (icode == kEMFSCObrems  || 
+         icode == kEMFSCOmoller || 
+         icode == kEMFSCObhabha || 
+         icode == kEMFSCOcompton ) 
+      {
+         px = fPint[0];
+         py = fPint[1];
+         pz = fPint[2];
+         e  = fPint[3];
+      } else if (icode == kKASKADdray  || 
+                icode == kKASKADbrems || 
+                icode == kKASKADpair) {
+         px = GENSTK.plr[0] * GENSTK.cxr[0];
+         py = GENSTK.plr[0] * GENSTK.cyr[0];
+         pz = GENSTK.plr[0] * GENSTK.czr[0];
+         e  = GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6];
+      } else {
+         Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) *  ParticleMassFPC(TRACKR.jtrack));
+         px = p*TRACKR.cxtrck;
+         py = p*TRACKR.cytrck;
+         pz = p*TRACKR.cztrck;
+         e  = TRACKR.etrack;
+      }
   }
   else
-    Warning("TrackMomentum","momentum not available");
+      Warning("TrackMomentum","momentum not available");
 }
 
 //______________________________________________________________________________ 
@@ -1507,11 +1638,12 @@ Double_t TFluka::TrackStep() const
 {
 // Return the length in centimeters of the current step
 // TRACKR.ctrack = total curved path
-  FlukaCallerCode_t caller = GetCaller();
-  if (caller == kBXEntering || caller == kBXExiting || 
-      caller == kENDRAW     || caller == kUSDRAW || 
-      caller == kUSTCKV     || caller == kMGResumedTrack)
-    return 0.0;
+    FlukaCallerCode_t caller = GetCaller();
+    if (caller == kBXEntering || caller == kBXExiting || 
+       caller == kENDRAW     || caller == kUSDRAW || 
+       caller == kUSTCKV     || caller == kMGResumedTrack ||
+       caller == kSODRAW)
+       return 0.0;
   else if (caller == kMGDRAW)
     return TRACKR.ctrack;
   else {
@@ -1531,8 +1663,10 @@ Double_t TFluka::TrackLength() const
     return TRACKR.cmtrck;
   else if (caller == kMGResumedTrack) 
     return TRACKR.spausr[8];
+  else if (caller == kSODRAW)
+      return 0.0;
   else {
-    Warning("TrackLength", "track length not available");
+    Warning("TrackLength", "track length not available for caller %5d \n", caller);
     return 0.0;
   } 
 }
@@ -1549,6 +1683,9 @@ Double_t TFluka::TrackTime() const
     return TRACKR.atrack;
   else if (caller == kMGResumedTrack)
     return TRACKR.spausr[3];
+  else if (caller == kSODRAW) {
+      return (FLKSTK.agestk[FLKSTK.npflka]);
+  }
   else {
     Warning("TrackTime", "track time not available");
     return 0.0;
@@ -1574,7 +1711,9 @@ Double_t TFluka::Edep() const
   FlukaCallerCode_t caller = GetCaller();
     
   if (caller == kBXExiting || caller == kBXEntering || 
-      caller == kUSDRAW    || caller == kMGResumedTrack) return 0.0;
+      caller == kUSDRAW    || caller == kMGResumedTrack ||
+      caller == kSODRAW) 
+      return 0.0;
   Double_t sum = 0;
   Int_t i = -1;
   
@@ -1635,9 +1774,12 @@ Int_t TFluka::TrackPid() const
 // Return the id of the particle transported
 // TRACKR.jtrack = identity number of the particle
   FlukaCallerCode_t caller = GetCaller();
-  if (caller != kEEDRAW) {
+  if (caller != kEEDRAW && caller != kSODRAW) {
      return PDGFromId( CorrectFlukaId() );
   }
+  else if (caller == kSODRAW) {
+      return PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
+  }
   else
     return -1000;
 }
@@ -1650,8 +1792,12 @@ Double_t TFluka::TrackCharge() const
 // TRACKR.jtrack = identity number of the particle
     
   FlukaCallerCode_t caller = GetCaller();
-  if (caller != kEEDRAW) 
-     return PAPROP.ichrge[CorrectFlukaId()+6];
+  if (caller != kEEDRAW && caller != kSODRAW) 
+     return PAPROP.ichrge[CorrectFlukaId() + 6];
+  else if (caller == kSODRAW) {
+      Int_t ifl =  PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
+      return PAPROP.ichrge[ifl + 6];
+  }
   else
     return -1000.0;
 }
@@ -1662,8 +1808,12 @@ Double_t TFluka::TrackMass() const
 // PAPROP.am = particle mass in GeV
 // TRACKR.jtrack = identity number of the particle
   FlukaCallerCode_t caller = GetCaller();
-  if (caller != kEEDRAW)
+  if (caller != kEEDRAW && caller != kSODRAW)
      return PAPROP.am[CorrectFlukaId()+6];
+  else if (caller == kSODRAW) {
+      Int_t ifl =  FLKSTK.iloflk[FLKSTK.npflka];
+      return PAPROP.am[ifl + 6];
+  }
   else
     return -1000.0;
 }
@@ -1672,11 +1822,34 @@ Double_t TFluka::TrackMass() const
 Double_t TFluka::Etot() const
 {
 // TRACKR.etrack = total energy of the particle
-  FlukaCallerCode_t caller = GetCaller();
-  if (caller != kEEDRAW)
-    return TRACKR.etrack;
-  else
-    return -1000.0;
+  FlukaCallerCode_t  caller = GetCaller();
+  FlukaProcessCode_t icode  = GetIcode();
+  if (caller != kEEDRAW && caller != kSODRAW && caller != kUSDRAW)
+  {
+      return TRACKR.etrack;
+  } else if (caller == kUSDRAW) {
+      if (icode == kEMFSCObrems  || 
+         icode == kEMFSCOmoller || 
+         icode == kEMFSCObhabha || 
+         icode == kEMFSCOcompton ) {
+         return  fPint[3];
+      }
+      else if (icode == kKASKADdray  || 
+              icode == kKASKADbrems || 
+              icode == kKASKADpair) {
+         return (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);      
+      }
+  }
+  else if (caller == kSODRAW) {
+      Int_t ist  = FLKSTK.npflka;
+      Double_t p = FLKSTK.pmoflk[ist];
+      Int_t ifl  = FLKSTK.iloflk[ist];
+      Double_t m = PAPROP.am[ifl + 6];
+      Double_t e = TMath::Sqrt(p * p + m * m);
+      return e;
+  }
+  
+  return -1000.0;
 }
 
 //
@@ -1917,39 +2090,69 @@ Int_t TFluka::StepProcesses(TArrayI &proc) const
   //
   // Return processes active in the current step
   //
-    FlukaProcessCode_t icode = GetIcode();
+    FlukaProcessCode_t icode   = GetIcode();
+    FlukaCallerCode_t  caller  = GetCaller();
+    
     proc.Set(1);
     TMCProcess iproc;
-    switch (icode) {
-    case kKASKADtimekill:
-    case kEMFSCOtimekill:
-    case kKASNEUtimekill:
-    case kKASHEAtimekill:
-    case kKASOPHtimekill:
-        iproc =  kPTOFlimit;
-        break;
-    case kKASKADstopping:
-    case kKASKADescape:
-    case kEMFSCOstopping1:
-    case kEMFSCOstopping2:
-    case kEMFSCOescape:
-    case kKASNEUstopping:
-    case kKASNEUescape:
-    case kKASHEAescape:
-    case kKASOPHescape:
-        iproc = kPStop;
-        break;
-    case kKASOPHabsorption:
-        iproc = kPLightAbsorption;
-        break;
-    case kKASOPHrefraction:
-        iproc = kPLightRefraction;
-    case kEMFSCOlocaldep : 
-        iproc = kPPhotoelectric;
-        break;
-    default:
-        iproc = ProdProcess(0);
+    if (caller == kBXEntering || caller == kBXExiting || caller == kEEDRAW || caller == kSODRAW) {
+       iproc = kPTransportation;
     }
+    else if (caller == kUSTCKV) {
+       iproc = kPCerenkov;
+    } else {
+       switch (icode) {
+       case kEMFSCO:
+           if (Edep() > 0.) {
+               iproc = kPEnergyLoss;
+           } else {
+               iproc = kPTransportation;
+           }
+           break;
+       case kKASKAD:
+           if (Edep() > 0.) {
+               iproc = kPEnergyLoss;
+           } else {
+               iproc = kPTransportation;
+           }
+           break;
+       case kKASHEA:
+       case kKASNEU:
+       case kKASOPH:
+       case kKASKADescape:
+       case kEMFSCOescape:
+       case kKASNEUescape:
+       case kKASHEAescape:
+       case kKASOPHescape:
+           iproc = kPTransportation;
+           break;
+       case kKASKADtimekill:
+       case kEMFSCOtimekill:
+       case kKASNEUtimekill:
+       case kKASHEAtimekill:
+       case kKASOPHtimekill:
+           iproc =  kPTOFlimit;
+           break;
+       case kKASKADstopping:
+       case kEMFSCOstopping1:
+       case kEMFSCOstopping2:
+       case kKASNEUstopping:
+           iproc = kPStop;
+           break;
+       case kKASOPHabsorption:
+           iproc = kPLightAbsorption;
+           break;
+       case kKASOPHrefraction:
+           iproc = kPLightRefraction;
+           break;
+       case kEMFSCOlocaldep : 
+           iproc = kPPhotoelectric;
+           break;
+       default:
+           iproc = ProdProcess(0);
+       }
+    }
+    
     proc[0] = iproc;
     return 1;
 }
@@ -1971,6 +2174,19 @@ const char* TFluka::VolName(Int_t id) const
    return fMCGeo->VolName(id);
 }
 
+Int_t TFluka::MediumId(const Text_t* mediumName) const
+{
+    //
+    // Return the unique medium id for medium with name mediumName
+    TList *medlist = gGeoManager->GetListOfMedia();
+    TGeoMedium* med = (TGeoMedium*) medlist->FindObject(mediumName);
+    if (med) {
+       return (med->GetId());
+    } else {
+       return (-1);
+    }
+}
+
 //______________________________________________________________________________ 
 Int_t TFluka::VolId(const Text_t* volName) const
 {
@@ -2021,7 +2237,7 @@ const char* TFluka::CurrentVolName() const
 //
 // Return the current volume name
 //
-  if (gGeoManager->IsOutside()) return 0;
+  if (gGeoManager->IsOutside()) return "OutOfWorld";
   return gGeoManager->GetCurrentVolume()->GetName();
 }
 
@@ -2337,8 +2553,7 @@ Double_t TFluka::GetPrimaryElectronKineticEnergy(Int_t i) const
     Double_t ekin = -1.;
     
     if (i >= 0 && i < ALLDLT.nalldl) {
-       Int_t j = ALLDLT.nalldl - 1 - i;
-        ekin =  ALLDLT.talldl[j];
+        ekin =  ALLDLT.talldl[i];
     } else {
         Warning("GetPrimaryElectronKineticEnergy",
                 "Primary electron index out of range %d %d \n",
@@ -2351,10 +2566,9 @@ void TFluka::GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Doubl
 {
     // Returns position  of primary electron i
         if (i >= 0 && i < ALLDLT.nalldl) {
-           Int_t j = ALLDLT.nalldl - 1 - i;
-           x = ALLDLT.xalldl[j];
-           y = ALLDLT.yalldl[j];
-           z = ALLDLT.zalldl[j];
+           x = ALLDLT.xalldl[i];
+           y = ALLDLT.yalldl[i];
+           z = ALLDLT.zalldl[i];
            return;
        } else {
            Warning("GetPrimaryElectronPosition",
@@ -2395,3 +2609,26 @@ void  TFluka::PrimaryIonisationStepping(Int_t nprim)
     // Reset the index
     SetCurrentPrimaryElectronIndex(-1);
 }
+
+//______________________________________________________________________
+Float_t* TFluka::CreateFloatArray(Double_t* array, Int_t size) const
+{
+// Converts Double_t* array to Float_t*,
+// !! The new array has to be deleted by user.
+// ---
+
+  Float_t* floatArray;
+  if (size>0) {
+    floatArray = new Float_t[size];
+    for (Int_t i=0; i<size; i++)
+      if (array[i] >= FLT_MAX ) 
+        floatArray[i] = FLT_MAX/100.;
+      else     
+        floatArray[i] = array[i];
+  }
+  else {
+    //floatArray = 0;
+    floatArray = new Float_t[1];
+  }
+  return floatArray;
+}