From 088f8e792ce3da25a5271d1d98a4d250d376e929 Mon Sep 17 00:00:00 2001
From: cholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Date: Fri, 9 Dec 2005 09:34:35 +0000
Subject: [PATCH] Fixed some coding style violations.

Added track lenght through medum to AliFMDHit.  Also added flag
fStop for stopped tracks (IsTrackStop || IsTrackDisappeared).

Updated DrawHits, DrawXsection.
Added GetMedia.C to track where tracks corresponding to hits come
from.  Useful for background studies.  Added DrawHitsDigits, to compare
the energy loss from hits to the ADC counts from digitisation.

Updated Config.C to properly deal with TFluka.
---
 FMD/AliFMD.cxx                      |  10 +-
 FMD/AliFMD.h                        |   7 +-
 FMD/AliFMD1.cxx                     |   3 +
 FMD/AliFMD1.h                       |  10 +-
 FMD/AliFMD2.cxx                     |   3 +
 FMD/AliFMD2.h                       |   6 +-
 FMD/AliFMD3.cxx                     |  15 ++-
 FMD/AliFMD3.h                       |   8 +-
 FMD/AliFMDBoolMap.cxx               |   6 +-
 FMD/AliFMDBoolMap.h                 |   2 +
 FMD/AliFMDDetector.cxx              |  59 ++++++++++-
 FMD/AliFMDDetector.h                |   2 +
 FMD/AliFMDDigit.h                   |  23 +++--
 FMD/AliFMDDigitizer.cxx             |  26 +++--
 FMD/AliFMDDigitizer.h               |   3 +
 FMD/AliFMDEdepMap.cxx               |  14 ++-
 FMD/AliFMDEdepMap.h                 |  14 ++-
 FMD/AliFMDGeometry.cxx              |   5 +-
 FMD/AliFMDGeometry.h                |   9 +-
 FMD/AliFMDGeometryBuilder.cxx       |   6 +-
 FMD/AliFMDGeometryBuilderSimple.cxx |   2 +-
 FMD/AliFMDHit.cxx                   |  27 +++--
 FMD/AliFMDHit.h                     |  12 ++-
 FMD/AliFMDMap.cxx                   |  19 ++--
 FMD/AliFMDMap.h                     |   3 +-
 FMD/AliFMDMult.cxx                  |   3 +-
 FMD/AliFMDMult.h                    |   8 +-
 FMD/AliFMDMultAlgorithm.cxx         |   6 +-
 FMD/AliFMDMultAlgorithm.h           |   4 +-
 FMD/AliFMDMultNaiive.cxx            |   3 +-
 FMD/AliFMDMultPoisson.cxx           |   2 +
 FMD/AliFMDMultRegion.h              |   8 +-
 FMD/AliFMDMultStrip.h               |   8 +-
 FMD/AliFMDParameters.cxx            |   2 +
 FMD/AliFMDParameters.h              |   3 +-
 FMD/AliFMDRawReader.h               |   4 +-
 FMD/AliFMDRawStream.h               |   1 +
 FMD/AliFMDRawWriter.h               |   3 +
 FMD/AliFMDRing.cxx                  |   7 ++
 FMD/AliFMDUShortMap.cxx             |   6 ++
 FMD/AliFMDv1.cxx                    |  36 ++++---
 FMD/AliFMDv1.h                      |   4 +-
 FMD/Config.C                        |  65 +++++++++++-
 FMD/scripts/.cvsignore              |   4 +
 FMD/scripts/DrawHits.C              |  70 +++++++------
 FMD/scripts/DrawHitsDigits.C        | 108 ++++++++++++++++++++
 FMD/scripts/DrawXsection.C          |  48 ++++++---
 FMD/scripts/GetMedia.C              | 153 ++++++++++++++++++++++++++++
 48 files changed, 693 insertions(+), 157 deletions(-)
 create mode 100644 FMD/scripts/DrawHitsDigits.C
 create mode 100644 FMD/scripts/GetMedia.C

diff --git a/FMD/AliFMD.cxx b/FMD/AliFMD.cxx
index 17fb122553c..50d9ac86001 100644
--- a/FMD/AliFMD.cxx
+++ b/FMD/AliFMD.cxx
@@ -236,6 +236,7 @@ AliFMD::~AliFMD ()
 AliFMD&
 AliFMD::operator=(const AliFMD& other)
 {
+  // Assignment operator
   AliDetector::operator=(other);
   fSDigits		= other.fSDigits; 
   fNsdigits		= other.fNsdigits;
@@ -830,7 +831,9 @@ AliFMD::AddHitByFields(Int_t    track,
 		       Float_t  pz,
 		       Float_t  edep,
 		       Int_t    pdg,
-		       Float_t  t)
+		       Float_t  t, 
+		       Float_t  l, 
+		       Bool_t   stop)
 {
   //
   // Add a hit to the list
@@ -851,6 +854,8 @@ AliFMD::AddHitByFields(Int_t    track,
   //    edep	  Energy deposited by track
   //    pdg	  Track's particle Id #
   //    t	  Time when the track hit 
+  //    l         Track length through the material. 
+  //    stop      Whether track was stopped or disappeared
   // 
   TClonesArray& a = *(HitsArray());
   // Search through the list of already registered hits, and see if we
@@ -877,7 +882,8 @@ AliFMD::AddHitByFields(Int_t    track,
   }
   // If hit wasn't already registered, do so know. 
   hit = new (a[fNhits]) AliFMDHit(fIshunt, track, detector, ring, sector, 
-				  strip, x, y, z, px, py, pz, edep, pdg, t);
+				  strip, x, y, z, px, py, pz, edep, pdg, t, 
+				  l, stop);
   fNhits++;
   return hit;
 }
diff --git a/FMD/AliFMD.h b/FMD/AliFMD.h
index c58f0d77f64..1b14ebdfddd 100644
--- a/FMD/AliFMD.h
+++ b/FMD/AliFMD.h
@@ -11,6 +11,9 @@
 //____________________________________________________________________
 //
 //  Manager class for the FMD - Base class.
+//  AliFMDv1, AliFMDv0, and AliFMDAlla 
+//  provides concrete implementations. 
+//  This class is sooooo crowded
 //
 #ifndef ALIDETECTOR_H 
 # include <AliDetector.h>
@@ -73,7 +76,9 @@ public:
 				       Float_t  pz=0,
 				       Float_t  edep=0,
 				       Int_t    pdg=0,
-				       Float_t  t=0);
+				       Float_t  t=0, 
+				       Float_t  len=0, 
+				       Bool_t   stopped=kFALSE);
   virtual        void   AddDigit(Int_t *digits, Int_t* notused=0);
   virtual        void   AddDigitByFields(UShort_t detector=0, 
 					 Char_t   ring='\0', 
diff --git a/FMD/AliFMD1.cxx b/FMD/AliFMD1.cxx
index bb2b038c622..98e2040b60f 100644
--- a/FMD/AliFMD1.cxx
+++ b/FMD/AliFMD1.cxx
@@ -20,6 +20,9 @@
 // Concrete implementation of AliFMDDetector 
 //
 // This implements the geometry for FMD1 
+// FMD1 has only one ring, of type `inner'.  
+// It is sitting at z=320.
+// It is the FMD ring with highest eta
 //
 #include "AliFMD1.h"		// ALIFMD1_H 
 #include "AliFMDRing.h"		// ALIFMDRING_H 
diff --git a/FMD/AliFMD1.h b/FMD/AliFMD1.h
index 29c295962cd..79a798cd37b 100644
--- a/FMD/AliFMD1.h
+++ b/FMD/AliFMD1.h
@@ -9,11 +9,11 @@
 class AliFMDRing;
 
 //__________________________________________________________________
-/** Geometry description and parameters of the FMD1
-    detector. 
-    
-    The FMD1 only has one ring.     
-*/
+// Geometry description and parameters of the FMD1
+// detector. 
+//
+// The FMD1 only has one ring.     
+//
 class AliFMD1 : public AliFMDDetector 
 {
 public:
diff --git a/FMD/AliFMD2.cxx b/FMD/AliFMD2.cxx
index 37ac44a271b..dcf5967c92f 100644
--- a/FMD/AliFMD2.cxx
+++ b/FMD/AliFMD2.cxx
@@ -20,6 +20,9 @@
 // Concrete implementation of AliFMDDetector 
 //
 // This implements the geometry for FMD2
+// The FMD2 has two ring, one of both types. 
+// FMD2 is mounted on the space-frame via 4 flanges
+// Support is not fleshed ot yet. 
 //
 #include "AliFMD2.h"		// ALIFMD2_H 
 #include "AliFMDRing.h"		// ALIFMDRING_H 
diff --git a/FMD/AliFMD2.h b/FMD/AliFMD2.h
index b0a0711a302..138c26b5855 100644
--- a/FMD/AliFMD2.h
+++ b/FMD/AliFMD2.h
@@ -8,9 +8,9 @@
 # include "AliFMDDetector.h"
 #endif
 
-/** Geometry description and parameters of the FMD2
-    detector. 
-*/
+// Geometry description and parameters of the FMD2
+// detector. 
+// This has two rings. 
 class AliFMD2 : public AliFMDDetector 
 {
 protected: 
diff --git a/FMD/AliFMD3.cxx b/FMD/AliFMD3.cxx
index 82fdf251fd9..83efcab8633 100644
--- a/FMD/AliFMD3.cxx
+++ b/FMD/AliFMD3.cxx
@@ -19,7 +19,10 @@
 //                                                                          
 // Concrete implementation of AliFMDDetector 
 //
-// This implements the geometry for FMD3
+// This implements the geometry for FMD3.
+// This has 2 rings.
+// The support of the FMD3 is a carbon-fibre cone, attached to the ITS
+// support via flanges.  The cone also supports the beam-pipe.
 //
 #include "AliFMD3.h"		// ALIFMD3_H 
 #include "AliLog.h"		// ALILOG_H
@@ -36,6 +39,7 @@ ClassImp(AliFMD3)
 AliFMD3::AliFMD3(AliFMDRing* inner, AliFMDRing* outer) 
   : AliFMDDetector(3, inner, outer)
 {
+  // Constructor. 
   SetInnerZ(-62.8);
   SetOuterZ(-75.2);
   SetNoseZ();
@@ -57,6 +61,7 @@ AliFMD3::AliFMD3(AliFMDRing* inner, AliFMDRing* outer)
 void
 AliFMD3::Init() 
 {
+  // Initialize 
   AliFMDDetector::Init();
   SetInnerHoneyHighR(GetOuterHoneyHighR());
   Double_t zdist   = fConeLength - fBackLength - fNoseLength;
@@ -76,16 +81,16 @@ AliFMD3::ConeR(Double_t z, Option_t* opt) const
 {
   // Calculate the cone radius at Z
   if (fAlpha < 0) {
-    Warning("ConeR", "alpha not set: %lf", fAlpha);
+    AliWarning(Form("alpha not set: %lf", fAlpha));
     return -1;
   }
   if (z > fNoseZ) {
-    Warning("ConeR", "z=%lf is before start of cone %lf", z, fNoseZ);
+    AliWarning(Form("z=%lf is before start of cone %lf", z, fNoseZ));
     return -1;
   }
   if (z < fOuterZ - fOuter->GetRingDepth() - fHoneycombThickness) {
-    Warning("ConeR", "z=%lf is after end of cone %lf", z, 
-	    fOuterZ - fOuter->GetRingDepth() - fHoneycombThickness);
+    AliWarning(Form("z=%lf is after end of cone %lf", z, 
+		    fOuterZ - fOuter->GetRingDepth() - fHoneycombThickness));
     return -1;
   }
   Double_t e = fBeamThickness / TMath::Cos(TMath::ATan(fAlpha));
diff --git a/FMD/AliFMD3.h b/FMD/AliFMD3.h
index 9a24496b0f6..e2292571950 100644
--- a/FMD/AliFMD3.h
+++ b/FMD/AliFMD3.h
@@ -8,10 +8,10 @@
 # include "AliFMDDetector.h"
 #endif
 
-/** Geometry description and parameters of the FMD3 detector.
-    
-    FMD3 has a fairly complicated support structure 
-*/
+// Geometry description and parameters of the FMD3 detector.
+// FMD3 has a fairly complicated support structure. 
+// The cone also supports the beam-pipe.
+//
 class AliFMD3 : public AliFMDDetector 
 {
 public: 
diff --git a/FMD/AliFMDBoolMap.cxx b/FMD/AliFMDBoolMap.cxx
index 918a7cbc70b..c3b413fd968 100644
--- a/FMD/AliFMDBoolMap.cxx
+++ b/FMD/AliFMDBoolMap.cxx
@@ -17,8 +17,10 @@
 /* $Id$ */
 //__________________________________________________________
 // 
-// Map of per strip Bool_t information
-// 
+// Map of Bool_t for each FMD strip
+// Used in calibration and the like classes.
+// Used amoung other things for dead-channel map
+// Can also be used for other stuff too
 // Created Mon Nov  8 12:51:51 2004 by Christian Holm Christensen
 // 
 #include "AliFMDBoolMap.h"	//ALIFMDBOOLMAP_H
diff --git a/FMD/AliFMDBoolMap.h b/FMD/AliFMDBoolMap.h
index a81892b86b6..61b00e6b878 100644
--- a/FMD/AliFMDBoolMap.h
+++ b/FMD/AliFMDBoolMap.h
@@ -10,6 +10,8 @@
 //__________________________________________________________
 // 
 // Map of Bool_t for each FMD strip
+// Used in calibration and the like classes.
+// Used amoung other things for dead-channel map
 // 
 #ifndef ALIFMDMAP_H
 # include <AliFMDMap.h>
diff --git a/FMD/AliFMDDetector.cxx b/FMD/AliFMDDetector.cxx
index e28f52c10c1..aaf59abe34d 100644
--- a/FMD/AliFMDDetector.cxx
+++ b/FMD/AliFMDDetector.cxx
@@ -27,7 +27,6 @@
 //
 #include "AliFMDDetector.h"	// ALIFMDSUBDETECTOR_H
 #include "AliFMDRing.h"		// ALIFMDRING_H
-#include <AliLog.h>		// ALILOG_H
 
 //====================================================================
 ClassImp(AliFMDDetector)
@@ -42,6 +41,12 @@ AliFMDDetector::AliFMDDetector(Int_t id, AliFMDRing* inner, AliFMDRing* outer)
     fInner(inner),
     fOuter(outer)
 {
+  // Constructor
+  // 
+  //   ID         Id of detector (1,2, or 3)
+  //   INNER      Inner ring geometry 
+  //   OUTER      Outer ring geometry (if any)
+  // 
   SetHoneycombThickness();
   SetAlThickness();
   SetInnerHoneyLowR(0);
@@ -52,10 +57,50 @@ AliFMDDetector::AliFMDDetector(Int_t id, AliFMDRing* inner, AliFMDRing* outer)
   SetOuterHoneyHighR(0);
 }
 
+//____________________________________________________________________
+AliFMDDetector::AliFMDDetector(const AliFMDDetector& other)
+  : TNamed(other), 
+    fId(other.fId),
+    fInner(other.fInner),
+    fOuter(other.fOuter)
+{
+  // Copy constructor 
+  SetHoneycombThickness(other.GetHoneycombThickness());
+  SetAlThickness(other.GetAlThickness());
+  SetInnerHoneyLowR(other.GetInnerHoneyLowR());
+  SetInnerHoneyHighR(other.GetInnerHoneyHighR());
+  SetInnerZ(other.GetInnerZ());
+  SetOuterZ(other.GetOuterZ());
+  SetOuterHoneyLowR(other.GetOuterHoneyLowR());
+  SetOuterHoneyHighR(other.GetOuterHoneyHighR());
+}
+
+//____________________________________________________________________
+AliFMDDetector&
+AliFMDDetector::operator=(const AliFMDDetector& other)
+{
+  // Assignment operator
+  SetName(other.GetName());
+  SetTitle(other.GetTitle());
+  fId    = other.fId;
+  fInner = other.fInner;
+  fOuter = other.fOuter;
+  SetHoneycombThickness(other.GetHoneycombThickness());
+  SetAlThickness(other.GetAlThickness());
+  SetInnerHoneyLowR(other.GetInnerHoneyLowR());
+  SetInnerHoneyHighR(other.GetInnerHoneyHighR());
+  SetInnerZ(other.GetInnerZ());
+  SetOuterZ(other.GetOuterZ());
+  SetOuterHoneyLowR(other.GetOuterHoneyLowR());
+  SetOuterHoneyHighR(other.GetOuterHoneyHighR());
+  return *this;
+}
+
 //____________________________________________________________________
 void
 AliFMDDetector::Init()
 {
+  // Initialize. 
   if (fInner) {
     SetInnerHoneyLowR(fInner->GetLowR() + 1.);
     SetInnerHoneyHighR(fInner->GetHighR() + 1.);
@@ -71,6 +116,10 @@ AliFMDDetector::Init()
 AliFMDRing*
 AliFMDDetector::GetRing(Char_t id) const
 {
+  // Get the specified ring 
+  // 
+  //   ID      Id of ring ('I' or 'O')
+  // 
   switch (id) {
   case 'i':
   case 'I': return GetInner();
@@ -84,6 +133,10 @@ AliFMDDetector::GetRing(Char_t id) const
 Double_t
 AliFMDDetector::GetRingZ(Char_t id) const
 {
+  // Get the z-coordinate specified ring 
+  // 
+  //   ID      Id of ring ('I' or 'O')
+  // 
   switch (id) {
   case 'i':
   case 'I': return GetInnerZ();
@@ -101,6 +154,8 @@ AliFMDDetector::Detector2XYZ(Char_t ring,
 			     Double_t& y, 
 			     Double_t& z) const
 {
+  // Translate detector coordinates (this,ring,sector,strip) into
+  // (x,y,z) coordinates (in global reference frame)
   AliFMDRing* r = GetRing(ring);
   if (!r) return;
   z = GetRingZ(ring);
@@ -116,6 +171,8 @@ AliFMDDetector::XYZ2Detector(Double_t  x,
 			     UShort_t& sector,
 			     UShort_t& strip) const
 {
+  // Translate (x,y,z) coordinates (in global reference frame) into 
+  // detector coordinates (this,ring,sector,strip).
   AliFMDRing* rng = 0;
   ring = -1;
   for (int j = 0; j < 2; j++) {
diff --git a/FMD/AliFMDDetector.h b/FMD/AliFMDDetector.h
index d4442ce7b81..3b4f7e8cd36 100644
--- a/FMD/AliFMDDetector.h
+++ b/FMD/AliFMDDetector.h
@@ -23,6 +23,8 @@ class AliFMDDetector : public TNamed
 {
 public:
   AliFMDDetector(Int_t id, AliFMDRing* inner, AliFMDRing* outer);
+  AliFMDDetector(const AliFMDDetector& other);
+  AliFMDDetector& operator=(const AliFMDDetector& other);
   virtual ~AliFMDDetector() {}
   /** Initialize the geometry */
   virtual void Init();
diff --git a/FMD/AliFMDDigit.h b/FMD/AliFMDDigit.h
index 42cde9138da..ba2ef9df17c 100644
--- a/FMD/AliFMDDigit.h
+++ b/FMD/AliFMDDigit.h
@@ -2,7 +2,10 @@
 #define ALIFMDDIGIT_H
 //___________________________________________________________________
 //
-//  Digits classes for the FMD                
+//  Digits classes for the FMD
+//  AliFMDBaseDigit - base class 
+//  AliFMDDigit     - Normal (smeared) digit             
+//  AliFMDSDigit    - Summable (non-smeared) digit             
 //
 #ifndef ROOT_TObject
 # include <TObject.h>
@@ -35,10 +38,6 @@ protected:
 //____________________________________________________________________
 class AliFMDDigit : public AliFMDBaseDigit
 {
-protected:
-  UShort_t fCount1;     // Digital signal 
-  Short_t  fCount2;     // Digital signal (-1 if not used)
-  Short_t  fCount3;     // Digital signal (-1 if not used)
 public:
   AliFMDDigit();
   AliFMDDigit(UShort_t detector, 
@@ -54,6 +53,10 @@ public:
   Short_t  Count3()                const { return fCount3;   }
   UShort_t Counts()                const;
   void     Print(Option_t* opt="") const;
+protected:
+  UShort_t fCount1;     // Digital signal 
+  Short_t  fCount2;     // Digital signal (-1 if not used)
+  Short_t  fCount3;     // Digital signal (-1 if not used)
   ClassDef(AliFMDDigit,1)     // Normal FMD digit
 };
 
@@ -68,11 +71,6 @@ AliFMDDigit::Counts() const
 //____________________________________________________________________
 class AliFMDSDigit : public AliFMDBaseDigit
 {
-protected:
-  Float_t  fEdep;       // Energy deposited 
-  UShort_t fCount1;     // Digital signal 
-  Short_t  fCount2;     // Digital signal (-1 if not used)
-  Short_t  fCount3;     // Digital signal (-1 if not used)
 public:
   AliFMDSDigit();
   AliFMDSDigit(UShort_t detector, 
@@ -90,6 +88,11 @@ public:
   Float_t  Edep()                  const { return fEdep;     }
   UShort_t Counts()                const;
   void     Print(Option_t* opt="") const;
+protected:
+  Float_t  fEdep;       // Energy deposited 
+  UShort_t fCount1;     // Digital signal 
+  Short_t  fCount2;     // Digital signal (-1 if not used)
+  Short_t  fCount3;     // Digital signal (-1 if not used)
   ClassDef(AliFMDSDigit,1)     // Summable FMD digit
 };
   
diff --git a/FMD/AliFMDDigitizer.cxx b/FMD/AliFMDDigitizer.cxx
index 46851ec3f86..e2ea608e104 100644
--- a/FMD/AliFMDDigitizer.cxx
+++ b/FMD/AliFMDDigitizer.cxx
@@ -440,25 +440,24 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t   edep,
   // 
   //                  = E + (l - E) * ext(-B * t)
   // 
-  const Float_t mipEnergy = 1.664 * siThickness * siDensity;
-  const Float_t convf     = (1 / mipEnergy * Float_t(fAltroChannelSize) 
-			     / fVA1MipRange);
-  UShort_t ped            = MakePedestal();
+  Float_t  mipEnergy  = 1.664 * siThickness * siDensity;
+  Float_t  convF      = (1/mipEnergy*Float_t(fAltroChannelSize)/fVA1MipRange);
+  UShort_t ped       = MakePedestal();
 
   // In case we don't oversample, just return the end value. 
   if (fSampleRate == 1) {
-    counts[0] = UShort_t(TMath::Min(edep * convf + ped, 
+    counts[0] = UShort_t(TMath::Min(edep * convF + ped, 
 				    Float_t(fAltroChannelSize)));
     return;
   }
-
+  
   // Create a pedestal 
   Int_t   n = fSampleRate;
-  Float_t B = fShapingTime;
+  Float_t b = fShapingTime;
   for (Ssiz_t i = 0; i < n;  i++) {
     Float_t t = Float_t(i) / n;
-    Float_t s = edep + (last - edep) * TMath::Exp(-B * t);
-    counts[i] = UShort_t(TMath::Min(s * convf + ped, 
+    Float_t s = edep + (last - edep) * TMath::Exp(-b * t);
+    counts[i] = UShort_t(TMath::Min(s * convF + ped, 
 				    Float_t(fAltroChannelSize))); 
   }
 }
@@ -558,6 +557,7 @@ AliFMDDigitizer::Exec(Option_t*)
 UShort_t
 AliFMDDigitizer::MakePedestal() const 
 {
+  // Make a pedestal 
   return UShort_t(TMath::Max(gRandom->Gaus(fPedestal, fPedestalWidth), 0.));
 }
 
@@ -573,6 +573,7 @@ AliFMDDigitizer::AddDigit(AliFMD*  fmd,
 			  Short_t  count2, 
 			  Short_t  count3) const
 {
+  // Add a digit
   fmd->AddDigitByFields(detector, ring, sector, strip, count1, count2, count3);
 }
 
@@ -582,14 +583,15 @@ AliFMDDigitizer::CheckDigit(Float_t         /* edep */,
 			    UShort_t        nhits,
 			    const TArrayI&  counts) 
 {
+  // Check that digit is consistent
   Int_t integral = counts[0];
   if (counts[1] >= 0) integral += counts[1];
   if (counts[2] >= 0) integral += counts[2];
   integral -= Int_t(fPedestal + 2 * fPedestalWidth);
   if (integral < 0) integral = 0;
   
-  Float_t convf = Float_t(fVA1MipRange) / fAltroChannelSize;
-  Float_t mips  = integral * convf;
+  Float_t convF = Float_t(fVA1MipRange) / fAltroChannelSize;
+  Float_t mips  = integral * convF;
   if (mips > Float_t(nhits) + .5 || mips < Float_t(nhits) - .5) 
     Warning("CheckDigit", "Digit -> %4.2f MIPS != %d +/- .5 hits", 
 	    mips, nhits);
@@ -628,6 +630,7 @@ AliFMDSDigitizer::AliFMDSDigitizer(const Char_t* headerFile,
 //____________________________________________________________________
 AliFMDSDigitizer::~AliFMDSDigitizer() 
 {
+  // Destructor
   AliLoader* loader = fRunLoader->GetLoader("FMDLoader");
   loader->CleanSDigitizer();
 }
@@ -689,6 +692,7 @@ AliFMDSDigitizer::AddDigit(AliFMD*  fmd,
 			   Short_t  count2, 
 			   Short_t  count3) const
 {
+  // Add a summable digit
   fmd->AddSDigitByFields(detector, ring, sector, strip, edep, 
 			 count1, count2, count3); 
 }
diff --git a/FMD/AliFMDDigitizer.h b/FMD/AliFMDDigitizer.h
index f80ce3ce8d8..81050669727 100644
--- a/FMD/AliFMDDigitizer.h
+++ b/FMD/AliFMDDigitizer.h
@@ -77,6 +77,9 @@ protected:
   UShort_t      fSampleRate;       // Times the ALTRO samples pre-amp.
   Float_t       fShapingTime;      // Shaping profile parameter
   
+  AliFMDBaseDigitizer(const AliFMDBaseDigitizer& o) 
+    : AliDigitizer(o) {}
+  AliFMDBaseDigitizer& operator=(const AliFMDBaseDigitizer&) { return *this; }
   ClassDef(AliFMDBaseDigitizer,0) // Base class for FMD digitizers
 };
 
diff --git a/FMD/AliFMDEdepMap.cxx b/FMD/AliFMDEdepMap.cxx
index e6238b6b52a..dea63beb2f8 100644
--- a/FMD/AliFMDEdepMap.cxx
+++ b/FMD/AliFMDEdepMap.cxx
@@ -33,6 +33,7 @@ AliFMDEdepMap::AliFMDEdepMap(const AliFMDEdepMap& other)
 	      other.fMaxStrips), 
     fData(0)
 {
+  // Copy constructor 
   fData = new AliFMDEdepHitPair[fMaxDetectors * fMaxRings * 
 				fMaxSectors * fMaxStrips];
   for (size_t i = 0; i < fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips;
@@ -64,6 +65,7 @@ AliFMDEdepMap::AliFMDEdepMap(size_t maxDet,
 AliFMDEdepMap&
 AliFMDEdepMap::operator=(const AliFMDEdepMap& other) 
 {
+  // Assignment operator
   fMaxDetectors = other.fMaxDetectors;
   fMaxRings     = other.fMaxRings;
   fMaxSectors   = other.fMaxSectors;
@@ -76,12 +78,22 @@ AliFMDEdepMap::operator=(const AliFMDEdepMap& other)
   return *this;
 }
 
+//____________________________________________________________________
+void
+AliFMDEdepMap::Reset() 
+{
+  // Reset to zero
+  for (size_t i = 0; i < fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips;
+       i++) { fData[i].fEdep = 0; fData[i].fN = 0; };
+}
+
 //____________________________________________________________________
 void
 AliFMDEdepMap::Reset(const AliFMDEdepHitPair& val) 
 {
+  // Reset to val
   for (size_t i = 0; i < fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips;
-       i++) fData[i] = val;
+       i++) { fData[i].fEdep = val.fEdep; fData[i].fN = val.fN; };
 }
 
 //____________________________________________________________________
diff --git a/FMD/AliFMDEdepMap.h b/FMD/AliFMDEdepMap.h
index 608b75f2cba..e9831a1f469 100644
--- a/FMD/AliFMDEdepMap.h
+++ b/FMD/AliFMDEdepMap.h
@@ -10,16 +10,21 @@
 #endif 
 //____________________________________________________________________
 //
-// Cache of Energy deposited, hit information perr strip
+// Cache of Energy deposited, hit information per strip.
+// Contains a pair of fEdep and fN
+// fEdep is the summed energy deposition, and fN is the number of hits
 //
 
 //____________________________________________________________________
 class AliFMDEdepHitPair 
 {
 public:
-  Float_t  fEdep;
-  UShort_t fN;
+  Float_t  fEdep; // summed energy deposition
+  UShort_t fN;    // Number of hits
   AliFMDEdepHitPair() : fEdep(0), fN(0) {}
+  AliFMDEdepHitPair& operator=(const AliFMDEdepHitPair& o) 
+  { fEdep = o.fEdep; fN    = o.fN; return *this; }
+  AliFMDEdepHitPair(const AliFMDEdepHitPair& o) : fEdep(o.fEdep), fN(o.fN) {}
 };
 
 //____________________________________________________________________
@@ -33,7 +38,8 @@ public:
 		size_t maxStr = kMaxStrips);
   virtual ~AliFMDEdepMap() { delete [] fData; }
   AliFMDEdepMap& operator=(const AliFMDEdepMap& other);
-  virtual void Reset(const AliFMDEdepHitPair& val=AliFMDEdepHitPair());
+  virtual void Reset();
+  virtual void Reset(const AliFMDEdepHitPair& val);
   virtual AliFMDEdepHitPair& operator()(UShort_t detector, 
 				     Char_t   ring, 
 				     UShort_t sector, 
diff --git a/FMD/AliFMDGeometry.cxx b/FMD/AliFMDGeometry.cxx
index 8743159288f..83031a45638 100644
--- a/FMD/AliFMDGeometry.cxx
+++ b/FMD/AliFMDGeometry.cxx
@@ -193,7 +193,10 @@ void
 AliFMDGeometry::SetActive(Int_t* active, Int_t n) 
 {
   fActive.Set(n);
-  for (Int_t i = 0; i < n; i++) fActive[i] = active[i];
+  for (Int_t i = 0; i < n; i++) { 
+    AliDebug(1, Form("Active vol id # %d: %d", i, active[i]));
+    fActive[i] = active[i];
+  }
 }
 
 //____________________________________________________________________
diff --git a/FMD/AliFMDGeometry.h b/FMD/AliFMDGeometry.h
index bee945a8272..210d9f44e6d 100644
--- a/FMD/AliFMDGeometry.h
+++ b/FMD/AliFMDGeometry.h
@@ -1,9 +1,14 @@
 #ifndef ALIFMDGEOMETRY_H
 #define ALIFMDGEOMETRY_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights
+ * reserved. 
+ *
+ * Latest changes by Christian Holm Christensen <cholm@nbi.dk>
+ *
+ * See cxx source for full Copyright notice                               
+ */
 //____________________________________________________________________ 
 //  
-// $Id$ 
-//
 // Forward Multiplicity Detector based on Silicon wafers. 
 //
 // This class is a singleton that handles the geometry parameters of
diff --git a/FMD/AliFMDGeometryBuilder.cxx b/FMD/AliFMDGeometryBuilder.cxx
index 4d3b2cd8746..93e2bbed3b7 100644
--- a/FMD/AliFMDGeometryBuilder.cxx
+++ b/FMD/AliFMDGeometryBuilder.cxx
@@ -131,7 +131,7 @@ AliFMDGeometryBuilder::AliFMDGeometryBuilder()
     fPlastic(0)
 {
   // Default constructor
-  fActiveId.Set(4);
+  fActiveId.Set(2);
 }
 
 //____________________________________________________________________
@@ -153,7 +153,7 @@ AliFMDGeometryBuilder::AliFMDGeometryBuilder(Bool_t detailed)
   //      fmd		Pointer to AliFMD object 
   //      detailed      Whether to make a detailed simulation or not 
   // 
-  fActiveId.Set(4);
+  fActiveId.Set(2);
 }
 
 
@@ -246,7 +246,7 @@ AliFMDGeometryBuilder::RingGeometry(AliFMDRing* r)
   case 'i':
   case 'I': fActiveId[0] = sid; break;
   case 'o':
-  case 'O': fActiveId[2] = sid; break;
+  case 'O': fActiveId[1] = sid; break;
   }
 
   // Shape of Printed circuit Board 
diff --git a/FMD/AliFMDGeometryBuilderSimple.cxx b/FMD/AliFMDGeometryBuilderSimple.cxx
index 2dff552a3bb..dc654af08d0 100644
--- a/FMD/AliFMDGeometryBuilderSimple.cxx
+++ b/FMD/AliFMDGeometryBuilderSimple.cxx
@@ -177,7 +177,7 @@ AliFMDGeometryBuilderSimple::RingGeometry(AliFMDRing* r)
   case 'i':
   case 'I': fActiveId[0] = sid; break;
   case 'o':
-  case 'O': fActiveId[2] = sid; break;
+  case 'O': fActiveId[1] = sid; break;
   }
 
   // Shape of Printed circuit Board 
diff --git a/FMD/AliFMDHit.cxx b/FMD/AliFMDHit.cxx
index 41bfc7b1c06..980cd2c9ea7 100644
--- a/FMD/AliFMDHit.cxx
+++ b/FMD/AliFMDHit.cxx
@@ -18,8 +18,15 @@
 //____________________________________________________________________
 //
 //  Hits in the FMD 
-//
-// Latest changes by Christian Holm Christensen
+//  Contains information on:
+//	Position of hit
+//	Momentum of track
+//	PID of track
+//	Energy loss of track
+//	Track #
+//	Track path length
+//	Track stopping status. 
+//  Latest changes by Christian Holm Christensen
 //
 #include "AliFMDHit.h"		// ALIFMDHIT_H
 #include "AliLog.h"		// ALILOG_H
@@ -46,8 +53,11 @@ AliFMDHit::AliFMDHit()
     fPz(0),
     fPdg(0),
     fEdep(0), 
-    fTime(0)
+    fTime(0), 
+    fLength(0), 
+    fStop(0)
 {
+  // Default CTOR
   fX = fY = fZ = 0;
 }
   
@@ -67,7 +77,9 @@ AliFMDHit::AliFMDHit(Int_t    shunt,
 		     Float_t  pz,
 		     Float_t  edep,
 		     Int_t    pdg,
-		     Float_t  t)
+		     Float_t  t, 
+		     Float_t  l, 
+		     Bool_t   stop)
   : AliHit(shunt, track),
     fDetector(detector), 
     fRing(ring), 
@@ -78,7 +90,9 @@ AliFMDHit::AliFMDHit(Int_t    shunt,
     fPz(pz),
     fPdg(pdg),
     fEdep(edep), 
-    fTime(t)
+    fTime(t), 
+    fLength(l), 
+    fStop(stop)
 {
   // Normal FMD hit ctor
   // 
@@ -152,7 +166,8 @@ AliFMDHit::Print(Option_t* option) const
     cout << "\tPDG:\t" << fPdg << " " << (pdg ? pdg->GetName() : "?") << "\n"
 	 << "\tP:\t(" << fPx << "," << fPy << "," << fPz << ") "<<P() << "\n" 
 	 << "\tX:\t" << fX << "," << fY << "," << fZ << "\n" 
-	 << "\tTrack #:\t" << fTrack << std::endl;
+	 << "\tTrack #:\t" << fTrack << "\tLength:\t" 
+	 << fLength << "cm\t" << (IsStop() ? "stopped" : "") << std::endl;
   }
 }
 
diff --git a/FMD/AliFMDHit.h b/FMD/AliFMDHit.h
index 8cddebe48a7..6e654047006 100644
--- a/FMD/AliFMDHit.h
+++ b/FMD/AliFMDHit.h
@@ -35,7 +35,9 @@ public:
 	    Float_t  pz=0,
 	    Float_t  edep=0,
 	    Int_t    pdg=0,
-	    Float_t  t=0);
+	    Float_t  t=0, 
+	    Float_t  l=0, 
+	    Bool_t   stop=kFALSE);
   virtual ~AliFMDHit() {}
 
   UShort_t Detector()	const { return fDetector; }
@@ -51,6 +53,8 @@ public:
   Float_t  Q()          const;
   Int_t    Pdg()        const { return fPdg;      }
   Float_t  Time()       const { return fTime;     }
+  Float_t  Length()     const { return fLength;   }
+  Bool_t   IsStop()     const { return fStop;     }
   void     Print(Option_t* opt="") const;
 
   void     SetEdep(Float_t edep) { fEdep = edep; }
@@ -65,8 +69,10 @@ protected:
   Int_t    fPdg;       // Particles PDG code 
   Float_t  fEdep;      // Energy deposition
   Float_t  fTime;      // Particle's time of flight
-
-  ClassDef(AliFMDHit,1)  //Hits for detector FMD
+  Float_t  fLength;    // Track length through material. 
+  Bool_t   fStop;      // Whether track was stopped. 
+  
+  ClassDef(AliFMDHit,2)  //Hits for detector FMD
 };
 #endif
 //____________________________________________________________________
diff --git a/FMD/AliFMDMap.cxx b/FMD/AliFMDMap.cxx
index 9b9149b4741..245088e231b 100644
--- a/FMD/AliFMDMap.cxx
+++ b/FMD/AliFMDMap.cxx
@@ -17,9 +17,14 @@
 
 //____________________________________________________________________
 //                                                                          
-//
+// Base class for caches of per-strip information.
+// This is used to index a strip. 
+// Data stored depends on derived class. 
+// This class provides some common infra-structure.
+// Derived classes sould define Reset, and operator(). 
 //
 #include "AliFMDMap.h"		// ALIFMDMAP_H
+#include "AliLog.h"
 
 //____________________________________________________________________
 ClassImp(AliFMDMap)
@@ -78,12 +83,12 @@ AliFMDMap::CalcIndex(size_t det, Char_t ring, size_t sec, size_t str) const
   Int_t idx = CheckIndex(det, ring, sec, str);
   if (idx < 0) {
     size_t ringi = (ring == 'I' ||  ring == 'i' ? 0 : 1);
-    Fatal("CalcIndex", "Index (%d,'%c',%d,%d) out of bounds, "
-	  "in particular the %s index ", 
-	  det, ring, sec, str, 
-	  (det >= fMaxDetectors ? "Detector" : 
-	   (ringi >= fMaxRings ? "Ring" : 
-	    (sec >= fMaxSectors ? "Sector" : "Strip"))));
+    AliFatal(Form("CalcIndex", "Index (%d,'%c',%d,%d) out of bounds, "
+		  "in particular the %s index ", 
+		  det, ring, sec, str, 
+		  (det >= fMaxDetectors ? "Detector" : 
+		   (ringi >= fMaxRings ? "Ring" : 
+		    (sec >= fMaxSectors ? "Sector" : "Strip")))));
     return 0;
   }
   return size_t(idx);
diff --git a/FMD/AliFMDMap.h b/FMD/AliFMDMap.h
index 4777ff908c6..2dba2c5529f 100644
--- a/FMD/AliFMDMap.h
+++ b/FMD/AliFMDMap.h
@@ -11,8 +11,9 @@
 //____________________________________________________________________
 //
 // Base class for caches of per-strip information.
+// This is used to index a strip. 
+// Data stored depends on derived class. 
 //
-
 class AliFMDMap : public TObject 
 {
 public:
diff --git a/FMD/AliFMDMult.cxx b/FMD/AliFMDMult.cxx
index 54f5e0edaf3..8a27f91598d 100644
--- a/FMD/AliFMDMult.cxx
+++ b/FMD/AliFMDMult.cxx
@@ -18,7 +18,8 @@
 //____________________________________________________________________
 //
 // Base class for reconstructed charged particle multiplicty in the
-// FMD.  The class contains the field fMethod which is a flag set by
+// FMD.  
+// The class contains the field fMethod which is a flag set by
 // the AliFMDMultAlgorithm that created the object. The flag tells us
 // which algorithm was used to create the data stored in the object. 
 //
diff --git a/FMD/AliFMDMult.h b/FMD/AliFMDMult.h
index ec4353741bc..a086e09afcc 100644
--- a/FMD/AliFMDMult.h
+++ b/FMD/AliFMDMult.h
@@ -1,10 +1,10 @@
 #ifndef ALIFMDMULT_H
 #define ALIFMDMULT_H
 
-/* Reconstracted Particles Class: has number of reconstructed
- * particles in sectors from NumOfMinSector to NumberOfMaxSector()
- * rings from NumOfMinRing to NumOfMaxRing for each FMDvolume
- */
+// Reconstracted Particles Class: has number of reconstructed
+// particles in sectors from NumOfMinSector to NumberOfMaxSector()
+// rings from NumOfMinRing to NumOfMaxRing for each FMDvolume 
+//
 #ifndef ROOT_TObject
 # include <TObject.h>
 #endif
diff --git a/FMD/AliFMDMultAlgorithm.cxx b/FMD/AliFMDMultAlgorithm.cxx
index 1e53582d486..7c367cff042 100644
--- a/FMD/AliFMDMultAlgorithm.cxx
+++ b/FMD/AliFMDMultAlgorithm.cxx
@@ -86,11 +86,14 @@ AliFMDMultAlgorithm::AliFMDMultAlgorithm(const char* name, const char* title)
     fTreeR(0), 
     fMult(0), 
     fFMD(0)
-{}
+{
+  // Default CTOR
+}
 
 //____________________________________________________________________
 AliFMDMultAlgorithm::~AliFMDMultAlgorithm()
 {
+  // DTOR
   if (fMult) {
     fMult->Delete();
     delete fMult;
@@ -102,6 +105,7 @@ AliFMDMultAlgorithm::~AliFMDMultAlgorithm()
 void
 AliFMDMultAlgorithm::PreEvent(TTree* treeR, Float_t /* ipZ */) 
 {
+  // Executed before each event.
   if (fMult) fMult->Clear();
   fNMult = 0;
   fTreeR = treeR;
diff --git a/FMD/AliFMDMultAlgorithm.h b/FMD/AliFMDMultAlgorithm.h
index 8719e91a8d7..8b074ddb300 100644
--- a/FMD/AliFMDMultAlgorithm.h
+++ b/FMD/AliFMDMultAlgorithm.h
@@ -8,10 +8,10 @@
  * See cxx source for full Copyright notice                               
  */
 /* $Id$ */
+
 //____________________________________________________________________
 // 
-// Base class for FMD reconstruction algorithms. 
-//
+// Base class for FMD reconstruction algorithms.
 // Derived classes will implement various ways of reconstructing the
 // charge particle multiplicity in the FMD.  
 // 
diff --git a/FMD/AliFMDMultNaiive.cxx b/FMD/AliFMDMultNaiive.cxx
index 412ffa2c9ee..1138fa3e89e 100644
--- a/FMD/AliFMDMultNaiive.cxx
+++ b/FMD/AliFMDMultNaiive.cxx
@@ -22,7 +22,8 @@
 // [See also the AliFMDReconstructor class]
 // 
 // This class reconstructs the multiplicity based on the assumption
-// that all particles are minimum ionizing. 
+// that all particles are minimum ionizing.   
+// Hence, the name `naiive'
 // 
 #include "AliFMD.h"			// ALIFMD_H
 #include "AliFMDMultNaiive.h"		// ALIFMDMULTNAIIVE_H
diff --git a/FMD/AliFMDMultPoisson.cxx b/FMD/AliFMDMultPoisson.cxx
index 18236d84b02..ef973d4ffd5 100644
--- a/FMD/AliFMDMultPoisson.cxx
+++ b/FMD/AliFMDMultPoisson.cxx
@@ -23,6 +23,7 @@
 // 
 // This class reconstructs the muliplicity in regions based on the
 // ratio of empty to full strips. 
+// Hence the name `poisson'
 //
 #include "AliFMD.h"			// ALIFMD_H
 #include "AliFMDGeometry.h"		// ALIFMDGEOMETRY_H
@@ -48,6 +49,7 @@ AliFMDMultPoisson::AliFMDMultPoisson()
     fDeltaPhi(0), 
     fThreshold(0)
 {
+  // CTOR
   SetDeltaEta();
   SetDeltaPhi();
   SetThreshold();
diff --git a/FMD/AliFMDMultRegion.h b/FMD/AliFMDMultRegion.h
index 46ad46c0767..af2a5090bb4 100644
--- a/FMD/AliFMDMultRegion.h
+++ b/FMD/AliFMDMultRegion.h
@@ -1,10 +1,10 @@
 #ifndef ALIFMDMULTREGION_H
 #define ALIFMDMULTREGION_H
 
-/* Reconstracted Particles Class: has number of reconstructed
- * particles in sectors from NumOfMinSector to NumberOfMaxSector()
- * rings from NumOfMinRing to NumOfMaxRing for each FMDvolume
- */
+// Reconstracted Particles Class: has number of reconstructed
+// particles in sectors from NumOfMinSector to NumberOfMaxSector()
+// rings from NumOfMinRing to NumOfMaxRing for each FMDvolume
+//
 #ifndef ALIFMDMULT_H
 # include "AliFMDMult.h"
 #endif
diff --git a/FMD/AliFMDMultStrip.h b/FMD/AliFMDMultStrip.h
index a2affc96195..d8d8b6b6824 100644
--- a/FMD/AliFMDMultStrip.h
+++ b/FMD/AliFMDMultStrip.h
@@ -1,10 +1,10 @@
 #ifndef ALIFMDMULTSTRIP_H
 #define ALIFMDMULTSTRIP_H
 
-/* Reconstracted Particles Class: has number of reconstructed
- * particles in sectors from NumOfMinSector to NumberOfMaxSector()
- * rings from NumOfMinRing to NumOfMaxRing for each FMDvolume
- */
+// Reconstracted Particles Class: has number of reconstructed
+// particles in sectors from NumOfMinSector to NumberOfMaxSector()
+// rings from NumOfMinRing to NumOfMaxRing for each FMDvolume
+//
 #ifndef ALIFMDMULT_H
 # include "AliFMDMult.h"
 #endif
diff --git a/FMD/AliFMDParameters.cxx b/FMD/AliFMDParameters.cxx
index 8ec4a78e00e..b152a3a7a10 100644
--- a/FMD/AliFMDParameters.cxx
+++ b/FMD/AliFMDParameters.cxx
@@ -21,6 +21,8 @@
 //
 // This class is a singleton that handles various parameters of
 // the FMD detectors.  
+// Eventually, this class will use the Conditions DB to get the
+// various parameters, which code can then request from here.
 //                                                       
 #include "AliFMDParameters.h"	// ALIFMDPARAMETERS_H
 #include "AliFMDGeometry.h"	// ALIFMDGEOMETRY_H
diff --git a/FMD/AliFMDParameters.h b/FMD/AliFMDParameters.h
index 062440c5dcf..496d5af12b2 100644
--- a/FMD/AliFMDParameters.h
+++ b/FMD/AliFMDParameters.h
@@ -12,6 +12,7 @@
 //
 //  Singleton class to handle various parameters (not geometry) of the
 //  FMD
+//  Should get ata fromm Conditions DB.
 //
 #ifndef ROOT_TNamed
 # include <TNamed.h>
@@ -49,7 +50,7 @@ public:
 protected:
   AliFMDParameters();
   virtual ~AliFMDParameters() {}
-  static AliFMDParameters* fgInstance;
+  static AliFMDParameters* fgInstance; // Static singleton instance
   
   const Float_t fSiDeDxMip;        // MIP dE/dx in Silicon
   UShort_t      fVA1MipRange;      // # MIPs the pre-amp can do    
diff --git a/FMD/AliFMDRawReader.h b/FMD/AliFMDRawReader.h
index 13ac0682497..01d95825dc5 100644
--- a/FMD/AliFMDRawReader.h
+++ b/FMD/AliFMDRawReader.h
@@ -11,6 +11,8 @@
 //____________________________________________________________________
 // 
 // Class to read ADC values from a AliRawReader object. 
+// Note, that it uses an ALTRO reader, which is wrong. 
+// Perhaps we need to implement it our selves
 // 
 #ifndef ROOT_TTask
 # include <TTask.h>
@@ -26,7 +28,7 @@ class AliFMDRawReader : public TTask
 {
 public:
   AliFMDRawReader(AliRawReader* reader, TTree* array);
-
+  virtual ~AliFMDRawReader() {}
   virtual void Exec(Option_t* option="");
 protected:
   TTree*        fTree;       //! Pointer to tree to read into 
diff --git a/FMD/AliFMDRawStream.h b/FMD/AliFMDRawStream.h
index ae891e3afa7..1fd692f4d16 100644
--- a/FMD/AliFMDRawStream.h
+++ b/FMD/AliFMDRawStream.h
@@ -24,6 +24,7 @@ class AliFMDRawStream : public AliAltroRawStream
 {
 public:
   AliFMDRawStream(AliRawReader* reader, UShort_t sampleRate=0);
+  virtual ~AliFMDRawStream() {}
 
   Short_t Sector()      const { return fRow; }
   Char_t  Ring()        const { return (fSector == 0 ? 'I' : 'O'); }
diff --git a/FMD/AliFMDRawWriter.h b/FMD/AliFMDRawWriter.h
index 94a2a0e83ed..0d142a93f44 100644
--- a/FMD/AliFMDRawWriter.h
+++ b/FMD/AliFMDRawWriter.h
@@ -11,6 +11,8 @@
 //____________________________________________________________________
 // 
 // Class to writer ADC values to a Raw File
+// Uses general ALTRO class - which is wrong
+// Should make it right!
 //
 #ifndef ROOT_TTask
 # include <TTask.h>
@@ -27,6 +29,7 @@ class AliFMDRawWriter : public TTask
 {
 public:
   AliFMDRawWriter(AliFMD* fmd);
+  virtual ~AliFMDRawWriter() {}
 
   virtual void Exec(Option_t* option="");
 protected:
diff --git a/FMD/AliFMDRing.cxx b/FMD/AliFMDRing.cxx
index c2b0ae88e37..18533c21ffc 100644
--- a/FMD/AliFMDRing.cxx
+++ b/FMD/AliFMDRing.cxx
@@ -44,6 +44,7 @@ AliFMDRing::AliFMDRing(Char_t id)
     fId(id), 
     fVerticies(0)
 {
+  // CTOR
   SetBondingWidth();
   SetWaferRadius();
   SetSiThickness();
@@ -74,6 +75,7 @@ AliFMDRing::AliFMDRing(Char_t id)
 void
 AliFMDRing::Init()
 {
+  // Initialize 
   Double_t tanTheta  = TMath::Tan(fTheta * TMath::Pi() / 180.);
   Double_t tanTheta2 = TMath::Power(tanTheta,2);
   Double_t r2        = TMath::Power(fWaferRadius,2);
@@ -112,6 +114,7 @@ AliFMDRing::Init()
 TVector2*
 AliFMDRing::GetVertex(Int_t i) const
 {
+  // Get the i'th vertex of polygon shape
   return static_cast<TVector2*>(fVerticies.At(i));
 }
 
@@ -123,6 +126,8 @@ AliFMDRing::Detector2XYZ(UShort_t sector,
 			 Double_t& y, 
 			 Double_t& z) const
 {
+  // Translate detector coordinates (this,sector,strip) to global
+  // coordinates (x,y,z)
   if (sector >= GetNSectors()) {
     Error("Detector2XYZ", "Invalid sector number %d (>=%d) in ring %c", 
 	  sector, GetNSectors(), fId);
@@ -149,6 +154,8 @@ AliFMDRing::XYZ2Detector(Double_t  x,
 			 UShort_t& sector,
 			 UShort_t& strip) const
 {
+  // Translate global coordinates (x,y,z) to detector coordinates
+  // (this,sector,strip)
   sector = strip = 0;
   Double_t r = TMath::Sqrt(x * x + y * y);
   Int_t str = Int_t((r - fMinR) / GetPitch());
diff --git a/FMD/AliFMDUShortMap.cxx b/FMD/AliFMDUShortMap.cxx
index 7404a47b351..a4be648921d 100644
--- a/FMD/AliFMDUShortMap.cxx
+++ b/FMD/AliFMDUShortMap.cxx
@@ -19,6 +19,9 @@
 //                                                                          
 // A map of per strip UShort_t information (for example ADC values,
 // number of hits and so on). 
+// Used for various calib information, and the like, 
+// as well as in reconstruction. 
+// Can be used elsewhere too.
 //
 #include "AliFMDUShortMap.h"		// ALIFMDUSHORTMAP_H
 
@@ -34,6 +37,7 @@ AliFMDUShortMap::AliFMDUShortMap(const AliFMDUShortMap& other)
 	      other.fMaxStrips), 
     fData(0)
 {
+  // CTOR
   fData = new UShort_t[fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips];
   for (size_t i = 0; i < fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips;
        i++) fData[i] = other.fData[i];
@@ -63,6 +67,7 @@ AliFMDUShortMap::AliFMDUShortMap(size_t maxDet,
 AliFMDUShortMap&
 AliFMDUShortMap::operator=(const AliFMDUShortMap& other) 
 {
+  // Assignment operator
   fMaxDetectors = other.fMaxDetectors;
   fMaxRings     = other.fMaxRings;
   fMaxSectors   = other.fMaxSectors;
@@ -78,6 +83,7 @@ AliFMDUShortMap::operator=(const AliFMDUShortMap& other)
 void
 AliFMDUShortMap::Reset(const UShort_t& val) 
 {
+  // Reset to val
   for (size_t i = 0; i < fMaxDetectors * fMaxRings * fMaxSectors * fMaxStrips;
        i++) fData[i] = val;
 }
diff --git a/FMD/AliFMDv1.cxx b/FMD/AliFMDv1.cxx
index e028469eccc..3155ec7c2ca 100644
--- a/FMD/AliFMDv1.cxx
+++ b/FMD/AliFMDv1.cxx
@@ -123,7 +123,11 @@ AliFMDv1::VMC2FMD(Int_t copy, TLorentzVector& v,
   detector = det;
 
   //Double_t  rz  = fmd->GetDetector(detector)->GetRingZ(ring);
-  Int_t     n   = fmd->GetDetector(detector)->GetRing(ring)->GetNSectors();
+  AliFMDDetector* gdet  = fmd->GetDetector(detector);
+  AliFMDRing*     gring = gdet->GetRing(ring);
+  if (!gring)
+    AliFatal(Form("Ring %c not found (%s)", ring, mc->CurrentVolPath()));
+  Int_t           n     = gring->GetNSectors();
 #if 0
   if (rz < 0) {
     Int_t s = ((n - sector + n / 2) % n) + 1;
@@ -133,7 +137,7 @@ AliFMDv1::VMC2FMD(Int_t copy, TLorentzVector& v,
   }
 #endif
   if (sector < 1 || sector > n) {
-    Warning("Step", "sector # %d out of range (0-%d)", sector-1, n-1);
+    AliWarning(Form("Step", "sector # %d out of range (0-%d)", sector-1, n-1));
     return kFALSE;
   }
   sector--;
@@ -239,17 +243,17 @@ AliFMDv1::StepManager()
     if (mc->IsTrackOut())         what.Append("out ");
     
     Int_t mother = gAlice->GetMCApp()->GetPrimary(trackno);
-    Warning("Step", "Track # %5d deposits a lot of energy\n" 
-	    "  Volume:    %s\n" 
-	    "  Momentum:  (%7.4f,%7.4f,%7.4f)\n"
-	    "  PDG:       %d (%s)\n" 
-	    "  Edep:      %-14.7f keV (mother %d)\n"
-	    "  p/m:       %-7.4f/%-7.4f = %-14.7f\n"
-	    "  Processes: %s\n"
-	    "  What:      %s\n",
-	    trackno, mc->CurrentVolPath(), p.X(), p.Y(), p.Z(),
-	    pdg, pname.Data(), edep, mother, p.P(), mass, 
-	    poverm, processes.Data(), what.Data());
+    AliWarning(Form("Step", "Track # %5d deposits a lot of energy\n" 
+		    "  Volume:    %s\n" 
+		    "  Momentum:  (%7.4f,%7.4f,%7.4f)\n"
+		    "  PDG:       %d (%s)\n" 
+		    "  Edep:      %-14.7f keV (mother %d)\n"
+		    "  p/m:       %-7.4f/%-7.4f = %-14.7f\n"
+		    "  Processes: %s\n"
+		    "  What:      %s\n",
+		    trackno, mc->CurrentVolPath(), p.X(), p.Y(), p.Z(),
+		    pdg, pname.Data(), edep, mother, p.P(), mass, 
+		    poverm, processes.Data(), what.Data()));
   }
   
   // Check that the track is actually within the active area 
@@ -286,11 +290,15 @@ AliFMDv1::StepManager()
 			"Accumulated Edep=%g (%f,%f,%f)", trackno, 
 			mc->CurrentVolPath(), edep, fCurrentDeltaE, 
 			v.X(), v.Y(), v.Z()));
+      TVector3 cur(v.Vect());
+      cur -= fCurrentV.Vect();
+      Double_t len = cur.Mag();
       AliFMDHit* h = 
 	AddHitByFields(trackno, detector, ring, sector, strip,
 		       fCurrentV.X(),  fCurrentV.Y(), fCurrentV.Z(),
 		       fCurrentP.X(),  fCurrentP.Y(), fCurrentP.Z(), 
-		       fCurrentDeltaE, fCurrentPdg,   fCurrentV.T());
+		       fCurrentDeltaE, fCurrentPdg,   fCurrentV.T(), 
+		       len, mc->IsTrackDisappeared()||mc->IsTrackStop());
       // Add a copy 
       if (isBad && fBad) { 
 	new ((*fBad)[fBad->GetEntries()]) AliFMDHit(*h);
diff --git a/FMD/AliFMDv1.h b/FMD/AliFMDv1.h
index 829d49d2c77..711c8e29da5 100644
--- a/FMD/AliFMDv1.h
+++ b/FMD/AliFMDv1.h
@@ -11,7 +11,9 @@
 //____________________________________________________________________
 //
 //  Manager class for the FMD - Detailed version. 
-//
+//  Implements the full geometry, 
+//  And does stepping
+// 
 #ifndef ALIFMD_H 
 # include "AliFMD.h"
 #endif
diff --git a/FMD/Config.C b/FMD/Config.C
index 0ac26284b51..c863044dc4e 100644
--- a/FMD/Config.C
+++ b/FMD/Config.C
@@ -113,6 +113,9 @@ enum EG_t {
   kMuonCocktailCent1Single,	//
   kMuonCocktailPer1Single,	//
   kMuonCocktailPer4Single,
+  kFlatFMD1, 
+  kFlatFMD2, 
+  kFlatFMD3,
   kEgMax
 };
 
@@ -171,7 +174,10 @@ const char* egName[kEgMax] = {
   "kMuonCocktailPer4HighPt",	//
   "kMuonCocktailCent1Single",	//
   "kMuonCocktailPer1Single",	//
-  "kMuonCocktailPer4Single"
+  "kMuonCocktailPer4Single",
+  "kFMD1Flat",
+  "kFMD2Flat",
+  "kFMD3Flat"
 };
 
 //____________________________________________________________________
@@ -316,10 +322,27 @@ Config()
 	((TFluka*)gMC)->SetGeneratePemf(kFALSE);
       else
 	((TFluka*)gMC)->SetGeneratePemf(kTRUE);
-      TString neuxsc(gSystem->Which(".", "neuxsc.bin"));
-      if (neusxc->IsNull()) {
-	gSystem->Link("$(FLUPRO)/neuxsc.bin", "neuxsc.bin");
+      TString flupro(gSystem->Getenv("FLUPRO"));
+      if (flupro.IsNull()) 
+	Fatal("Config.C", "Environment variable FLUPRO not set");
+#if 0
+      char* files[] = { "brems_fin.bin", "cohff.bin", "elasct.bin", 
+			"gxsect.bin", "nuclear.bin", "sigmapi.bin", 
+			0 };
+      char* file = files[0];
+      while (file) {
+	TString which(gSystem->Which(".", file));
+	if (which.IsNull()) {
+	  if (gSystem->Symlink(Form("%s/%s", flupro.Data(), file), file)!=0) 
+	    Fatal("Config.C", "Couldn't link $(FLUPRO)/%s -> .", file);
+	}
+	file++;
       }
+#endif
+      TString neuxsc(gSystem->Which(".", "neuxsc.bin"));
+      if (neuxsc.IsNull()) 
+	gSystem->Symlink(Form("%s/neuxsc_72.bin", flupro.Data()), 
+			 "neuxsc.bin"); 
       gSystem->CopyFile("$(FLUPRO)/random.dat", "old.seed", kTRUE);
     }
     break;
@@ -1551,6 +1574,40 @@ GeneratorFactory(EG_t eg, Rad_t rad, TString& comment)
       gGener=gener;
     }
     break;
+  case kFlatFMD1: 
+    {
+      comment = comment.Append(" Flat in FMD1 range");
+      AliGenBox* gener = AliGenBox(2000);
+      gener->SetPart(211);
+      gener->SetMomentumRange(3,4);
+      gener->SetPhiRange(0, 360);
+      gener->SetThetaRange(0.77, 3.08);
+      gGener = gener;
+    }
+    break;
+  case kFlatFMD2: 
+    {
+      comment = comment.Append(" Flat in FMD2 range");
+      AliGenBox* gener = AliGenBox(2000);
+      gener->SetPart(211);
+      gener->SetMomentumRange(3,4);
+      gener->SetPhiRange(0, 360);
+      gener->SetThetaRange(2.95, 20.42);
+      gGener = gener;
+    }
+    break;
+  case kFlatFMD3: 
+    {
+      comment = comment.Append(" Flat in FMD3 range");
+      AliGenBox* gener = AliGenBox(2000);
+      gener->SetPart(211);
+      gener->SetMomentumRange(3,4);
+      gener->SetPhiRange(0, 360);
+      gener->SetThetaRange(155.97, 176.73);
+      gGener = gener;
+    }
+    break;
+    
   default: break;
   }
   return gGener;
diff --git a/FMD/scripts/.cvsignore b/FMD/scripts/.cvsignore
index 3367afdbbf9..a6c30b1891c 100644
--- a/FMD/scripts/.cvsignore
+++ b/FMD/scripts/.cvsignore
@@ -1 +1,5 @@
 old
+*.so
+*.root
+peg
+
diff --git a/FMD/scripts/DrawHits.C b/FMD/scripts/DrawHits.C
index 2d2052ddded..eb82a1cc56a 100644
--- a/FMD/scripts/DrawHits.C
+++ b/FMD/scripts/DrawHits.C
@@ -21,32 +21,31 @@ class DrawHits : public AliFMDInputHits
 {
 private:
   TH2D* fElossVsPMQ; // Histogram 
+  Int_t fPdg;
 public:
-  DrawHits() 
-  { 
-    Double_t emax = 1e4;
-    Double_t emin = 1e-5;
-    Int_t    n    = 901;
-    TArrayF tkine(n);
-    Float_t dp   = 1/TMath::Log10(emax/emin)/10;
-    Float_t pmin = TMath::Log10(emin);
-    tkine[0]     = emin;
-    for (Int_t i=1; i < tkine.fN; i++) {
-      Float_t el = pmin + i * dp;
-      tkine[i]   = TMath::Power(10, el);
-    }
-    Double_t dmin = .1;
-    Double_t dmax = 10;
-    Int_t    nd = 1000;
-    TArrayF eloss(nd+1);
-    eloss[0] = dmin;
-    for (Int_t i = 1; i < eloss.fN; i++){
-      eloss[i] = dmin + i * (dmax-dmin)/(eloss.fN-1);
+  TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
+  {
+    TArrayF bins(n+1);
+    Float_t dp   = n / TMath::Log10(max / min);
+    Float_t pmin = TMath::Log10(min);
+    bins[0]      = min;
+    for (Int_t i = 1; i < n+1; i++) {
+      Float_t p = pmin + i / dp;
+      bins[i]   = TMath::Power(10, p);
     }
-    fElossVsPMQ = new TH2D("bad", "#Delta E vs. p/(mq^{2})", 
+    return bins;
+  }
+  DrawHits(Int_t m=1000, Double_t emin=1, Double_t emax=1000, 
+	   Int_t n=900, Double_t tmin=1e-2, Double_t tmax=1e3, 
+	   Int_t pdg=211) 
+    : fPdg(pdg)
+  { 
+    TArrayF tkine(MakeLogScale(n, tmin, tmax));
+    TArrayF eloss(MakeLogScale(m, emin, emax));
+    fElossVsPMQ = new TH2D("bad", "#Delta E/#Delta x vs. p/(mq^{2})", 
 			   tkine.fN-1, tkine.fArray, eloss.fN-1, eloss.fArray);
-    fElossVsPMQ->SetXTitle("p/(mq^{2})");
-    fElossVsPMQ->SetYTitle("#Delta E [MeV]");
+    fElossVsPMQ->SetXTitle(Form("p%s", (fPdg == 0 ? "/(mq^{2})" : " [GeV]")));
+    fElossVsPMQ->SetYTitle("#Delta E/#Delta x [MeV/cm]");
   }
   Bool_t ProcessHit(AliFMDHit* hit, TParticle*) 
   {
@@ -54,19 +53,30 @@ public:
       std::cout << "No hit" << std::endl;
       return kFALSE;
     }
-    if (TMath::Abs(hit->Pdg()) != 211) return kTRUE;
-    // Float_t pmq = 0;
-    // Float_t q = hit->Q() / 3.;
-    // if (hit->M() != 0 && hit->Q() != 0) 
-    //  pmq = hit->P() / (hit->M()*q*q);
-    Float_t pmq = hit->P();
-    fElossVsPMQ->Fill(pmq, hit->Edep());
+
+    if (hit->IsStop()) return kTRUE;
+    Float_t x = hit->P();
+    Float_t y = hit->Edep()/hit->Length();
+    if (fPdg != 0) {
+      if (TMath::Abs(hit->Pdg()) != fPdg) return kTRUE;
+    }
+    else {
+      Float_t q = hit->Q() / 3.;
+      if (hit->M() == 0 || q == 0) return kTRUE;
+      x /= hit->M();
+      y /= q * q;
+    }
+    fElossVsPMQ->Fill(x, y);
     return kTRUE;
   }
   Bool_t Finish()
   {
     gStyle->SetPalette(1);
     gStyle->SetOptTitle(0);
+    gStyle->SetCanvasColor(0);
+    gStyle->SetCanvasBorderSize(0);
+    gStyle->SetPadColor(0);
+    gStyle->SetPadBorderSize(0);
     fElossVsPMQ->SetStats(kFALSE);
     fElossVsPMQ->Draw("COLZ");
     return kTRUE;
diff --git a/FMD/scripts/DrawHitsDigits.C b/FMD/scripts/DrawHitsDigits.C
new file mode 100644
index 00000000000..3e43c80f0e7
--- /dev/null
+++ b/FMD/scripts/DrawHitsDigits.C
@@ -0,0 +1,108 @@
+//
+// $Id$
+//
+// Script that contains a class to draw eloss from hits, versus ADC
+// counts from digits, using the AliFMDInputHits class in the util library. 
+//
+// It draws the energy loss versus the p/(mq^2).  It can be overlayed
+// with the Bethe-Bloc curve to show how the simulation behaves
+// relative to the expected. 
+//
+// Use the script `Compile.C' to compile this class using ACLic. 
+//
+#include <TH2D.h>
+#include <AliFMDHit.h>
+#include <AliFMDDigit.h>
+#include <AliFMDInput.h>
+#include <AliFMDEdepMap.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TArrayF.h>
+
+class DrawHitsDigits : public AliFMDInputHits
+{
+private:
+  TH2D* fElossVsAdc; // Histogram 
+  AliFMDEdepMap fMap;
+public:
+  TArrayF MakeLogScale(Int_t n, Double_t min, Double_t max) 
+  {
+    TArrayF bins(n+1);
+    Float_t dp   = n / TMath::Log10(max / min);
+    Float_t pmin = TMath::Log10(min);
+    bins[0]      = min;
+    for (Int_t i = 1; i < n+1; i++) {
+      Float_t p = pmin + i / dp;
+      bins[i]   = TMath::Power(10, p);
+    }
+    return bins;
+  }
+  DrawHitsDigits(Int_t n=900, Double_t emin=1e-3, Double_t emax=10, 
+		 Int_t m=1100, Double_t amin=-0.5, Double_t amax=1099.5) 
+  { 
+    AddLoad(kDigits);
+    TArrayF eloss(MakeLogScale(n, emin, emax));
+    TArrayF adcs(m+1);
+    adcs[0] = amin;
+    for (Int_t i = 1; i < m+1; i++) adcs[i] = adcs[i-1] + (amax-amin)/m;
+    fElossVsAdc = new TH2D("bad", "#Delta E vs. ADC", 
+			   eloss.fN-1, eloss.fArray, adcs.fN-1, adcs.fArray);
+    fElossVsAdc->SetXTitle("#Delta E/#Delta x [MeV/cm]");
+    fElossVsAdc->SetYTitle("ADC value");
+  }
+  Bool_t Begin(Int_t ev) 
+  {
+    fMap.Reset();
+    return AliFMDInputHits::Begin(ev);
+  }
+  Bool_t ProcessHit(AliFMDHit* hit, TParticle*) 
+  {
+    // Cache the energy loss 
+    if (!hit) return kFALSE;
+    UShort_t det = hit->Detector();
+    Char_t   rng = hit->Ring();
+    UShort_t sec = hit->Sector();
+    UShort_t str = hit->Strip();
+    fMap(det, rng, sec, str).fEdep += hit->Edep();
+    fMap(det, rng, sec, str).fN++;
+    return kTRUE;
+  }
+  Bool_t Event() 
+  {
+    if (!AliFMDInputHits::Event()) return kFALSE;
+    Int_t nEv = fTreeD->GetEntries();
+    for (Int_t i = 0; i < nEv; i++) {
+      Int_t digitRead  = fTreeD->GetEntry(i);
+      if (digitRead <= 0) continue;
+      Int_t nDigit = fArrayD->GetEntries();
+      if (nDigit <= 0) continue;
+      for (Int_t j = 0; j < nDigit; j++) {
+	AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
+	if (!digit) continue;
+	UShort_t det = digit->Detector();
+	Char_t   rng = digit->Ring();
+	UShort_t sec = digit->Sector();
+	UShort_t str = digit->Strip();
+	fElossVsAdc->Fill(fMap(det, rng, sec, str).fEdep, digit->Counts());
+      }    
+    }
+    return kTRUE;
+  }
+  Bool_t Finish()
+  {
+    gStyle->SetPalette(1);
+    gStyle->SetOptTitle(0);
+    gStyle->SetCanvasColor(0);
+    gStyle->SetCanvasBorderSize(0);
+    gStyle->SetPadColor(0);
+    gStyle->SetPadBorderSize(0);
+    fElossVsAdc->SetStats(kFALSE);
+    fElossVsAdc->Draw("COLZ");
+    return kTRUE;
+  }
+};
+
+//____________________________________________________________________
+//
+// EOF
+//
diff --git a/FMD/scripts/DrawXsection.C b/FMD/scripts/DrawXsection.C
index 36cdab1dd23..c4388f70cf7 100644
--- a/FMD/scripts/DrawXsection.C
+++ b/FMD/scripts/DrawXsection.C
@@ -1,5 +1,9 @@
+//
+// Script to draw a X-section, LOSS, or range made with MakeXsection
+//
 void
-DrawXsection(const char* filename="xsec.root", 
+DrawXsection(Bool_t scale=kFALSE, 
+	     const char* filename="xsec.root", 
 	     const char* var="LOSS", 
 	     const char* medName="FMD_Si$", 
 	     Double_t thick=.03,
@@ -19,26 +23,38 @@ DrawXsection(const char* filename="xsec.root",
   vb->SetAddress(&value);
   Int_t n = tree->GetEntries();
 
-  TDatabasePDG* pdgDb = TDatabasePDG::Instance();
-  TParticlePDG* pdgP  = pdgDb->GetParticle(pdgName);
-  if (!pdgP) {
-    std::cerr << "Couldn't find particle " << pdgName << std::endl;
-    return;
+  Float_t xscale = 1;
+  Float_t yscale = 1;
+  if (scale) {
+    TDatabasePDG* pdgDb = TDatabasePDG::Instance();
+    TParticlePDG* pdgP  = pdgDb->GetParticle(pdgName);
+    if (!pdgP) {
+      std::cerr << "Couldn't find particle " << pdgName << std::endl;
+      return;
+    }
+    Double_t m = pdgP->Mass();
+    Double_t q = pdgP->Charge() / 3;
+    if (m == 0 || q == 0) {
+      std::cerr  << "Mass is 0" << std::endl;
+      return;
+    }
+    xscale = 1 / m;
+    yscale = 1 / (q * q);
   }
-  // Double_t m = pdgP->Mass();
-  // Double_t q = pdgP->Charge() / 3;
-  // std::cout << q << "\t" << m << std::endl;
-  // if (m == 0) {
-  ///  std::cerr  << "Mass is 0" << std::endl;
-  // return;
-  // }
   
-  TGraph* graph = new TGraph(n);
+  TGraphErrors* graph = new TGraphErrors(n);
   for (Int_t i = 0; i < n; i++) {
     tree->GetEntry(i);
-    graph->SetPoint(i, tkine, value*thick); // /(m*q*q)
+    Double_t x = tkine*xscale;
+    Double_t y = value*yscale;
+    graph->SetPoint(i, x, y); 
+    // 5 sigma
+    graph->SetPointError(i, 0, 5 * .1 * y);
   }
-  graph->Draw("LP same");
+  graph->SetLineWidth(2);
+  graph->SetFillStyle(3001);
+  graph->SetFillColor(6);
+  graph->Draw("L3 same");
 }
 
 //____________________________________________________________________
diff --git a/FMD/scripts/GetMedia.C b/FMD/scripts/GetMedia.C
new file mode 100644
index 00000000000..6186081bee8
--- /dev/null
+++ b/FMD/scripts/GetMedia.C
@@ -0,0 +1,153 @@
+//
+// $Id$
+//
+// Script that contains a class to get the media where a certain track
+// was created.   This is used for background studies. 
+//
+// Use the script `Compile.C' to compile this class using ACLic. 
+//
+#include <TH2D.h>
+#include <AliFMDHit.h>
+#include <AliFMDInput.h>
+#include <iostream>
+#include <TStyle.h>
+#include <TArrayF.h>
+#include <TArrayI.h>
+#include <AliRun.h>
+#include <TNtuple.h>
+
+struct Media : public TNamed
+{
+  TArrayI*        fMeds;
+  TNtuple*        fCount;
+  Media(const char* name, const char* title) 
+    : TNamed(name, title), fMeds(0) 
+  {
+    fCount = new TNtuple(GetName(), "E/I:DET:RNG:SEC:STR:X/F:Y:Z:EDEP:PDG/I");
+  }
+  Media(const char* name) 
+    : TNamed(name, "Media information"), fMeds(0), fCount(0)
+  {
+    AliDetector* det = gAlice->GetDetector(name);
+    if (!det) {
+      Warning("Media", "Detector %s not found in gAlice", name);
+      return;
+    }
+    fMeds = det->GetIdtmed();
+    if (!fMeds) {
+      Warning("Media", "No mediums for detector %s found", name);
+      return;
+    }
+    fCount = new TNtuple(GetName(), "E/I:DET:RNG:SEC:STR:X/F:Y:Z:EDEP:PDG/I");
+  }
+};
+
+
+class GetMedia : public AliFMDInputHits
+{
+private:
+  TString    fModList;
+  TObjArray  fMedia;
+  Media*     fOther;
+  Media*     fAll;
+  Int_t      fEv;
+  TFile*     fOutput;
+public:
+  GetMedia(const char* modlist="FMD:ITS:BODY:ABSO:START:PIPE", 
+	   const char* output="media.root") 
+  { 
+    fOutput = TFile::Open(output, "RECREATE");
+    fOther  = new Media("other", "Unknown media"),
+    fAll    = new Media("All", "All media")
+    fEv     = 0;
+    AddLoad(kKinematics);
+  }
+  Media* FindMedia(Int_t med) 
+  {
+    TIter next(&fMedia);
+    Media* media = 0;
+    while ((media == static_cast<Media*>(next()))) {
+      if (!media->fMeds) continue;
+      TArrayI& ids = *(media->fMeds);
+      for (Int_t i = 0; i < ids.fN; i++) 
+	if (med == ids[i]) return media;
+    }
+    return 0;
+  }
+  Bool_t Init()
+  {
+    if (!gGeoManager) {
+      Error("GetMedia", "No geometry defined - make sure you have that");
+      return kFALSE;
+    }
+    if (!gAlice) {
+      Error("GetMedia", "gAlice not defined");
+      return kFALSE;
+    }
+    Int_t now = 0;
+    Int_t colon;
+    while ((colon = fModList.Index(":", now)) >= 0) {
+      fMedia.Add(new Media(fModList(now, colon-now)));
+      now = colon+1;
+    }
+    if (now < fModList.Length()) 
+      fMedia.Add(new Media(now, fModList.Length()-now));
+    if (fMedia.GetEntries() <= 0) return kFALSE;
+    return AliFMDInputHits::Init();
+  }  
+  Bool_t Begin(Int_t ev) 
+  {
+    fEv = ev;
+    return AliFMDInputHits::Begin(ev);
+  }
+  Bool_t ProcessHit(AliFMDHit* hit, TParticle* track) 
+  {
+    if (!hit || !track) {
+      std::cout << "No hit or track (hit: " << hit << ", track: " 
+		<< track << ")" << std::endl;
+      return kFALSE;
+    }
+    // Get production vertex 
+    Double_t vx = track->Vx();
+    Double_t vy = track->Vy();
+    Double_t vz = track->Vz();
+    
+    fAll->fCounts->Fill(fEv, hit->Detector(),Int_t(hit->Ring()),
+			hit->Sector(), hit->Strip(), 
+			hit->X(), hit->Y(), hit->Z(), hit->Edep(), 
+			hit->Pdg());
+    // Get node 
+    TGeoNode* prodNode = gGeoManager->FindNode(vx,vy,vz);
+    if (!prodNode) return kTRUE;
+    
+    // Get volume 
+    TGeoVolume* prodVol = prodNode->GetVolume();
+    if (!prodVol) return kTRUE;
+    
+    // Med medium 
+    TGeoMedium* prodMed = prodVol->GetMedium();
+    if (!prodMed) return kTRUE;
+    
+    Media* media = FindMedia(prodMed->GetUniqueID());
+    if (media) media = fOther; 
+    
+    // r = TMath::Sqrt(hit->X() * hit->X() + hit->Y() * hit->Y());
+    // if(r!=0) Float_t wgt=1./(2. * 3.1415*r);
+    media->fCounts->Fill(fEv, hit->Detector(),Int_t(hit->Ring()),
+			 hit->Sector(), hit->Strip(), 
+			 hit->X(), hit->Y(), hit->Z(), hit->Edep(), 
+			 hit->Pdg());
+    return kTRUE;
+  }
+  Bool_t Finish()
+  {
+    fOutput->Write();
+    fOutput->Close();
+    return kTRUE;
+  }
+};
+
+//____________________________________________________________________
+//
+// EOF
+//
-- 
2.43.0