Corrected indices. (E. Futo)
[u/mrichter/AliRoot.git] / TFluka / TFluka.cxx
index c3631df5a1e619d02a736a865f48dff078a9be1f..e8aa67a429d2b07c3818e5732d4b7689747657d4 100644 (file)
 
 /*
 $Log$
+Revision 1.9  2002/12/06 12:41:29  morsch
+Mess from last merge cleaned up.
+
+Revision 1.8  2002/12/06 12:28:44  morsch
+Region to media mapping corrected and improved.
+
 Revision 1.7  2002/12/06 12:21:32  morsch
 User stepping methods added (E. Futo)
 
@@ -499,592 +505,32 @@ Int_t TFluka::PDGFromId(Int_t id) const
   Int_t intfluka = GetFlukaIPTOKP(id);
   //MPKDHA() goes from internal to PDG
   return mpdgha(intfluka);
-  
 }
-<<<<<<< TFluka.cxx
-
-
 
 //_____________________________________________________________________________
 // methods for step management
 //____________________________________________________________________________ 
 //
-// dynamic properties
+// set methods
 //
-void TFluka::TrackPosition(TLorentzVector& position) const
-{
-// Return the current position in the master reference frame of the
-// track being transported
-// TRACKR.atrack = age of the particle
-// TRACKR.xtrack = x-position of the last point
-// TRACKR.ytrack = y-position of the last point
-// TRACKR.ztrack = z-position of the last point
-  position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
-  position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
-  position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
-  position.SetT(TRACKR.atrack);
-}
-
-void TFluka::TrackMomentum(TLorentzVector& momentum) const
-{
-// Return the direction and the momentum (GeV/c) of the track
-// currently being transported
-// TRACKR.ptrack = momentum of the particle (not always defined, if
-//               < 0 must be obtained from etrack) 
-// TRACKR.cx,y,ztrck = direction cosines of the current particle
-// TRACKR.etrack = total energy of the particle
-// TRACKR.jtrack = identity number of the particle
-// PAPROP.am[TRACKR.jtrack] = particle mass in gev
-  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 - PAPROP.am[TRACKR.jtrack]*PAPROP.am[TRACKR.jtrack]);
-    momentum.SetPx(p*TRACKR.cxtrck);
-    momentum.SetPy(p*TRACKR.cytrck);
-    momentum.SetPz(p*TRACKR.cztrck);
-    momentum.SetE(TRACKR.etrack);
-    return;
-  }
-}
-
-Double_t TFluka::TrackStep() const
-{
-// Return the length in centimeters of the current step
-// TRACKR.ctrack = total curved path
-    return TRACKR.ctrack;
-}
-
-Double_t TFluka::TrackLength() const
-{
-// It is wrong
-// should be the sum of all steps starting from the beginning of the track
-// for the time being returns only the length in centimeters of the current step
-    return TRACKR.ctrack;
-}
-
-Double_t TFluka::TrackTime() const
-{
-// Return the current time of flight of the track being transported
-// TRACKR.atrack = age of the particle
-  return TRACKR.atrack;
-}
-
-Double_t TFluka::Edep() const
-{
-// Energy deposition
-// if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
-// -->local energy deposition (the value and the point are not recorded in TRACKR)
-//    but in the variable "rull" of the procedure "endraw.cxx"
-// if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
-// -->no energy loss along the track
-// if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
-// -->energy loss distributed along the track
-// TRACKR.dtrack = energy deposition of the jth deposition even
-  if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
-    return fRull;
-  else {
-    Double_t sum = 0;
-    for ( Int_t j=0;j<TRACKR.mtrack;j++) {
-      sum +=TRACKR.dtrack[j];  
-    }
-    return sum;
-  }
-}
-
-Int_t TFluka::TrackPid() const
-{
-// Return the id of the particle transported
-// TRACKR.jtrack = identity number of the particle
-  return PDGFromId(TRACKR.jtrack);
-}
-
-Double_t TFluka::TrackCharge() const
-{
-// Return charge of the track currently transported
-// PAPROP.ichrge = electric charge of the particle
-  return PAPROP.ichrge[TRACKR.jtrack+6];
-}
-
-Double_t TFluka::TrackMass() const
-{
-// PAPROP.am = particle mass in GeV
-  return PAPROP.am[TRACKR.jtrack+6];
-}
-
-Double_t TFluka::Etot() const
-{
-// TRACKR.etrack = total energy of the particle
-  return TRACKR.etrack;
-}
-
-//
-// track status
-//
-Bool_t   TFluka::IsNewTrack() const
-{
-// ???????????????,
-// True if the track is not at the boundary of the current volume
-  return 0;
-}
-
-Bool_t   TFluka::IsTrackInside() const
-{
-// True if the track is not at the boundary of the current volume
-// In Fluka a step is always inside one kind of material
-// If the step would go behind the region of one material,
-// it will be shortened to reach only the boundary.
-// Therefore IsTrackInside() is always true.
-  return 1;
-}
-
-Bool_t   TFluka::IsTrackEntering() const
+void TFluka::SetMaxStep(Double_t)
 {
-// True if this is the first step of the track in the current volume
-// Boundary- (X) crossing
-// Icode = 19: boundary crossing - call from Kaskad
-// Icode = 29: boundary crossing - call from Emfsco
-// Icode = 39: boundary crossing - call from Kasneu
-// Icode = 49: boundary crossing - call from Kashea
-// Icode = 59: boundary crossing - call from Kasoph
-  if (fIcode == 19 ||
-      fIcode == 29 ||
-      fIcode == 39 ||
-      fIcode == 49 ||
-      fIcode == 59) return 1;
-  else return 0;
+// SetMaxStep is dummy procedure in TFluka !
+  cout << "SetMaxStep is dummy procedure in TFluka !" << endl;
 }
 
-Bool_t   TFluka::IsTrackExiting() const
+void TFluka::SetMaxNStep(Int_t)
 {
-// True if this is the last step of the track in the current volume
-// Boundary- (X) crossing
-// Icode = 19: boundary crossing - call from Kaskad
-// Icode = 29: boundary crossing - call from Emfsco
-// Icode = 39: boundary crossing - call from Kasneu
-// Icode = 49: boundary crossing - call from Kashea
-// Icode = 59: boundary crossing - call from Kasoph
-  if (fIcode == 19 ||
-      fIcode == 29 ||
-      fIcode == 39 ||
-      fIcode == 49 ||
-      fIcode == 59) return 1;
-  else return 0;
+// SetMaxNStep is dummy procedure in TFluka !
+  cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
 }
 
-Bool_t   TFluka::IsTrackOut() const
+void TFluka::SetUserDecay(Int_t)
 {
-// True if the track is out of the setup
-// means escape
-// Icode = 14: escape - call from Kaskad
-// Icode = 23: escape - call from Emfsco
-// Icode = 32: escape - call from Kasneu
-// Icode = 40: escape - call from Kashea
-// Icode = 51: escape - call from Kasoph
-  if (fIcode == 14 ||
-      fIcode == 23 ||
-      fIcode == 32 ||
-      fIcode == 40 ||
-      fIcode == 51) return 1;
-  else return 0;
-}
-
-Bool_t   TFluka::IsTrackDisappeared() const
-{
-// means all inelastic interactions and decays
-// Icode = 11: inelastic interaction recoil - call from Kaskad
-  if (fIcode == 11) return 1;
-  else return 0;
+// SetUserDecay is dummy procedure in TFluka !
+  cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
 }
 
-Bool_t   TFluka::IsTrackStop() const
-{
-// True if the track energy has fallen below the threshold
-// means stopped by signal or below energy threshold
-// Icode = 12: stopping particle       - call from Kaskad
-// Icode = 15: time kill               - call from Kaskad
-// Icode = 21: below threshold, iarg=1 - call from Emfsco
-// Icode = 22: below threshold, iarg=2 - call from Emfsco
-// Icode = 24: time kill               - call from Emfsco
-// Icode = 31: below threshold         - call from Kasneu
-// Icode = 33: time kill               - call from Kasneu
-// Icode = 41: time kill               - call from Kashea
-// Icode = 52: time kill               - call from Kasoph
-  if (fIcode == 12 ||
-      fIcode == 15 ||
-      fIcode == 21 ||
-      fIcode == 22 ||
-      fIcode == 24 ||
-      fIcode == 31 ||
-      fIcode == 33 ||
-      fIcode == 41 ||
-      fIcode == 52) return 1;
-  else return 0;
-}
-
-Bool_t   TFluka::IsTrackAlive() const
-{
-// means not disappeared or not out
-  if (IsTrackDisappeared() || IsTrackOut() ) return 0;
-  else return 1;
-}
-
-//
-// secondaries
-//
-
-//Int_t NSecondaries() const
-// Number of secondary particles generated in the current step
-//{
-//  return FINUC.np;
-//}
-
-//void     GetSecondary(Int_t isec, Int_t& particleId, TLorentzVector& position,
-//                   TLorentzVector& momentum)
-// 
-//{
-// will come from FINUC when called from USDRAW
-//}
-
-//TMCProcess ProdProcess(Int_t isec) const
-// Name of the process that has produced the secondary particles
-// in the current step
-//{
-// will come from FINUC when called from USDRAW
-//}
-
-//Int_t StepProcesses(TArrayI &proc) const
-// Return processes active in the current step
-//{
-//ck = total energy of the particl ???????????????? 
-//}
-
-
-void TFluka::FutoTest() const
-{
-  Int_t icode, mreg, newreg, medium;
-  Double_t rull, xsco, ysco, zsco;
-  icode  = GetIcode();
-  medium = -1;
-  
-  if (icode == 0)
-    cout << " icode=" << icode << endl;
-  else if (icode > 0 && icode <= 5) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    cout << " icode=" << icode
-        << " mreg=" << mreg
-        << " medium=" << medium
-        << endl;
-  }
-  else if (icode >= 10 && icode <= 15) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    rull = GetRull();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " rull=" << rull
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode >= 20 && icode <= 24) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    rull = GetRull();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-        << " mreg=" << mreg
-        << " medium=" << medium
-        << " rull=" << rull
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode >= 30 && icode <= 33) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    rull = GetRull();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " rull=" << rull
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode >= 40 && icode <= 41) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    rull = GetRull();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " rull=" << rull
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode >= 50 && icode <= 52) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    rull = GetRull();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " rull=" << rull
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode >= 100 && icode <= 105) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 208) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 210) {
-    mreg = GetMreg();
-   medium = GetMedium();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-        << " mreg=" << mreg
-        << " medium=" << medium 
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 212) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode >= 214 && icode <= 215) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 217) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 219) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 221) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 225) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 300) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 400) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 19) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    newreg = GetNewreg();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " newreg=" << newreg
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 29) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    newreg = GetNewreg();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " newreg=" << newreg
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 39) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    newreg = GetNewreg();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " newreg=" << newreg
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 49) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    newreg = GetNewreg();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " newreg=" << newreg
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-  else if (icode == 59) {
-    mreg = GetMreg();
-    medium = GetMedium();
-    newreg = GetNewreg();
-    xsco = GetXsco();
-    ysco = GetYsco();
-    zsco = GetZsco();
-    cout << " icode=" << icode
-         << " mreg=" << mreg
-        << " medium=" << medium
-        << " newreg=" << newreg
-        << " xsco=" << xsco
-        << " ysco=" << ysco
-        << " zsco=" << zsco << endl;
-  }
-
-
-} // end of FutoTest
-
-=======
-
-
-
-//_____________________________________________________________________________
-// methods for step management
-//____________________________________________________________________________ 
 //
 // dynamic properties
 //
@@ -1120,7 +566,7 @@ void TFluka::TrackMomentum(TLorentzVector& momentum) const
     return;
   }
   else {
-    Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack]*PAPROP.am[TRACKR.jtrack]);
+    Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
     momentum.SetPx(p*TRACKR.cxtrck);
     momentum.SetPy(p*TRACKR.cytrck);
     momentum.SetPz(p*TRACKR.cztrck);
@@ -1138,10 +584,17 @@ Double_t TFluka::TrackStep() const
 
 Double_t TFluka::TrackLength() const
 {
-// It is wrong
-// should be the sum of all steps starting from the beginning of the track
+// Still wrong !!!
+// This is the sum of substeps !!!
+// TRACKR.ctrack = total curved path of the current step
+// Sum of the substeps is identical to TRACKR.ctrack if the is no mag. field
+// The sum of all step length starting from the beginning of the track
 // for the time being returns only the length in centimeters of the current step
-    return TRACKR.ctrack;
+    Double_t sum = 0;
+    for ( Int_t j=0;j<TRACKR.ntrack;j++) {
+      sum +=TRACKR.ttrack[j];
+    }
+    return sum;
 }
 
 Double_t TFluka::TrackTime() const
@@ -1184,12 +637,14 @@ Double_t TFluka::TrackCharge() const
 {
 // Return charge of the track currently transported
 // PAPROP.ichrge = electric charge of the particle
+// TRACKR.jtrack = identity number of the particle
   return PAPROP.ichrge[TRACKR.jtrack+6];
 }
 
 Double_t TFluka::TrackMass() const
 {
 // PAPROP.am = particle mass in GeV
+// TRACKR.jtrack = identity number of the particle
   return PAPROP.am[TRACKR.jtrack+6];
 }
 
@@ -1323,65 +778,117 @@ Bool_t   TFluka::IsTrackAlive() const
 
 Int_t TFluka::NSecondaries() const
 // Number of secondary particles generated in the current step
-// fIcode = 100 = elastic interaction
-// fIcode = 101 = inelastic interaction
-// fIcode = 102 = particle decay
-// fIcode = 103 = delta ray
-// fIcode = 104 = pair production
-// fIcode = 105 = bremsstrahlung
+// FINUC.np = number of secondaries except light and heavy ions
+// FHEAVY.npheav = number of secondaries for light and heavy secondary ions
 {
-  if (fIcode >= 100 && fIcode <= 105)
-    return FINUC.np + FHEAVY.npheav;
-  else 
-    return -1;
+  return FINUC.np + FHEAVY.npheav;
 }
 
 void     TFluka::GetSecondary(Int_t isec, Int_t& particleId,
                TLorentzVector& position, TLorentzVector& momentum)
-// 
-// fIcode = 100 = elastic interaction
-// fIcode = 101 = inelastic interaction
-// fIcode = 102 = particle decay
-// fIcode = 103 = delta ray
-// fIcode = 104 = pair production
-// fIcode = 105 = bremsstrahlung
 {
-  if (fIcode >= 100 && fIcode <= 105) {
-    if (isec >= 0 && isec < FINUC.np) {
-      particleId = PDGFromId(FINUC.kpart[isec]);
-      position.SetX(fXsco);
-      position.SetY(fYsco);
-      position.SetZ(fZsco);
-      position.SetT(FINUC.agesec[isec]);
-      momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
-      momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
-      momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
-      momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
-    }
-    if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
-      Int_t jsec = isec - FINUC.np;
-      particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
-      position.SetX(fXsco);
-      position.SetY(fYsco);
-      position.SetZ(fZsco);
-      position.SetT(FHEAVY.agheav[jsec]);
-      momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
-      momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
-      momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
-      if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6) 
-        momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
-      else if (FHEAVY.tkheav[jsec] > 6)
-        momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
-    }
+  if (isec >= 0 && isec < FINUC.np) {
+         // more fine condition depending on icode
+         // icode = 100 ?
+         // icode = 101 OK
+         // icode = 102 OK
+         // icode = 103 ?
+         // icode = 104 ?
+         // icode = 105 ?
+         // icode = 208 ?
+         // icode = 210 ?
+         // icode = 212 ?
+         // icode = 214 OK
+         // icode = 215 OK
+         // icode = 219 ?
+         // icode = 221 OK
+         // icode = 225 ?
+         // icode = 300 ?
+         // icode = 400 ?
+         
+    particleId = PDGFromId(FINUC.kpart[isec]);
+    position.SetX(fXsco);
+    position.SetY(fYsco);
+    position.SetZ(fZsco);
+    position.SetT(TRACKR.atrack);
+//  position.SetT(TRACKR.atrack+FINUC.agesec[isec]); //not yet implem.
+    momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
+    momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
+    momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
+    momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
   }
+  if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
+    Int_t jsec = isec - FINUC.np;
+    particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
+    position.SetX(fXsco);
+    position.SetY(fYsco);
+    position.SetZ(fZsco);
+    position.SetT(TRACKR.atrack);
+//  position.SetT(TRACKR.atrack+FHEAVY.agheav[jsec]); //not yet implem.
+    momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
+    momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
+    momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
+    if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6) 
+      momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
+    else if (FHEAVY.tkheav[jsec] > 6)
+      momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
+    }
 }
 
-//TMCProcess ProdProcess(Int_t isec) const
+TMCProcess TFluka::ProdProcess(Int_t isec) const
 // Name of the process that has produced the secondary particles
 // in the current step
-//{
-// will come from FINUC when called from USDRAW
-//}
+{
+  const TMCProcess kIpNoProc = kPNoProcess;
+  const TMCProcess kIpPDecay = kPDecay;
+  const TMCProcess kIpPPair = kPPair;
+//const TMCProcess kIpPPairFromPhoton = kPPairFromPhoton;
+//const TMCProcess kIpPPairFromVirtualPhoton = kPPairFromVirtualPhoton;
+  const TMCProcess kIpPCompton = kPCompton;
+  const TMCProcess kIpPPhotoelectric = kPPhotoelectric;
+  const TMCProcess kIpPBrem = kPBrem;
+//const TMCProcess kIpPBremFromHeavy = kPBremFromHeavy;
+//const TMCProcess kIpPBremFromElectronOrPositron = kPBremFromElectronOrPositron;
+  const TMCProcess kIpPDeltaRay = kPDeltaRay;
+//const TMCProcess kIpPMoller = kPMoller;
+//const TMCProcess kIpPBhabha = kPBhabha;
+  const TMCProcess kIpPAnnihilation = kPAnnihilation;
+//const TMCProcess kIpPAnnihilInFlight = kPAnnihilInFlight;
+//const TMCProcess kIpPAnnihilAtRest = kPAnnihilAtRest;
+  const TMCProcess kIpPHadronic = kPHadronic;
+  const TMCProcess kIpPMuonNuclear = kPMuonNuclear;
+  const TMCProcess kIpPPhotoFission = kPPhotoFission;
+  const TMCProcess kIpPRayleigh = kPRayleigh;
+  const TMCProcess kIpPCerenkov = kPCerenkov;
+  const TMCProcess kIpPSynchrotron = kPSynchrotron;
+
+  Int_t mugamma = TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11;
+  if (fIcode == 102) return kIpPDecay;
+  else if (fIcode == 104 || fIcode == 217) return kIpPPair;
+//else if (fIcode == 104) return kIpPairFromPhoton;
+//else if (fIcode == 217) return kIpPPairFromVirtualPhoton;
+  else if (fIcode == 219) return kIpPCompton;
+  else if (fIcode == 221) return kIpPPhotoelectric;
+  else if (fIcode == 105 || fIcode == 208) return kIpPBrem;
+//else if (fIcode == 105) return kIpPBremFromHeavy;
+//else if (fIcode == 208) return kPBremFromElectronOrPositron;
+  else if (fIcode == 103 || fIcode == 400) return kIpPDeltaRay;
+  else if (fIcode == 210 || fIcode == 212) return kIpPDeltaRay;
+//else if (fIcode == 210) return kIpPMoller;
+//else if (fIcode == 212) return kIpPBhabha;
+  else if (fIcode == 214 || fIcode == 215) return kIpPAnnihilation;
+//else if (fIcode == 214) return kIpPAnnihilInFlight;
+//else if (fIcode == 215) return kIpPAnnihilAtRest;
+  else if (fIcode == 101) return kIpPHadronic;
+  else if (fIcode == 101) {
+    if (!mugamma) return kIpPHadronic;
+    else if (TRACKR.jtrack == 7) return kIpPPhotoFission;
+    else return kIpPMuonNuclear;
+  }
+  else if (fIcode == 225) return kIpPRayleigh;
+// Fluka codes 100, 300 and 400 still to be investigasted
+  else return kIpNoProc;
+}
 
 //Int_t StepProcesses(TArrayI &proc) const
 // Return processes active in the current step
@@ -1506,8 +1013,6 @@ void TFluka::FutoTest()
     cout << "TrackPid=" << TrackPid() << endl;
     cout << "NSecondaries=" << NSecondaries() << endl;
     for (Int_t isec=0; isec< NSecondaries(); isec++) {
-//void     TFluka::GetSecondary(Int_t isec, Int_t& particleId,
-//                 TLorentzVector& position, TLorentzVector& momentum)
       TFluka::GetSecondary(isec, particleId, position, momentum);
       cout << "TLorentzVector positionX=" << position.X()
            << "positionY=" << position.Y()