New sub-structure for fitting energy loss spectra with
authorcholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 29 Nov 2010 21:55:19 +0000 (21:55 +0000)
committercholm <cholm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 29 Nov 2010 21:55:19 +0000 (21:55 +0000)
landaus convolved with gaussians.  See also NIM B1,16

PWG2/FORWARD/analysis2/AliForwardUtil.cxx
PWG2/FORWARD/analysis2/AliForwardUtil.h

index bdf2227f2c080508428bbde6fbcf46387477f3d5..1dae159df5a8cc326e3a030c0b9c27cbb4a50930 100644 (file)
@@ -9,7 +9,233 @@
 #include <AliMultiplicity.h>
 #include <TH2D.h>
 #include <TH1I.h>
+#include <TF1.h>
+#include <TFitResult.h>
 #include <TMath.h>
+#include <TError.h>
+
+//====================================================================
+Int_t    AliForwardUtil::fgConvolutionSteps  = 100;
+Double_t AliForwardUtil::fgConvolutionNSigma = 5;
+namespace {
+  /** 
+   * The shift of the most probable value for the ROOT function TMath::Landau 
+   */
+  const Double_t  mpshift  = -0.22278298;
+  /** 
+   * Integration normalisation 
+   */
+  const Double_t  invSq2pi = 1. / TMath::Sqrt(2*TMath::Pi());
+
+  /** 
+   * Utility function to use in TF1 defintition 
+   */
+  Double_t landauGaus1(Double_t* xp, Double_t* pp) 
+  {
+    Double_t x        = xp[0];
+    Double_t constant = pp[0];
+    Double_t delta    = pp[1];
+    Double_t xi       = pp[2];
+    Double_t sigma    = pp[3];
+    Double_t sigma_n  = pp[4];
+
+    return constant * AliForwardUtil::LandauGaus(x, delta, xi, sigma, sigma_n);
+  }
+
+  /** 
+   * Utility function to use in TF1 defintition 
+   */
+  Double_t landauGausN(Double_t* xp, Double_t* pp) 
+  {
+    Double_t  x        = xp[0];
+    Double_t  constant = pp[0];
+    Double_t  delta    = pp[1];
+    Double_t  xi       = pp[2];
+    Double_t  sigma    = pp[3];
+    Double_t  sigma_n  = pp[4];
+    Int_t     n        = Int_t(pp[5]);
+    Double_t* a        = &(pp[6]);
+
+    return constant * AliForwardUtil::NLandauGaus(x, delta, xi, sigma, sigma_n,
+                                                 n, a);
+  }
+
+
+}
+//____________________________________________________________________
+Double_t 
+AliForwardUtil::Landau(Double_t x, Double_t delta, Double_t xi)
+{
+  return TMath::Landau(x, delta - xi * mpshift, xi);
+}
+//____________________________________________________________________
+Double_t 
+AliForwardUtil::LandauGaus(Double_t x, Double_t delta, Double_t xi,
+                          Double_t sigma, Double_t sigma_n)
+{
+  Double_t deltap = delta - xi * mpshift;
+  Double_t sigma2 = sigma_n*sigma_n + sigma*sigma;
+  Double_t sigma1 = TMath::Sqrt(sigma2);
+  Double_t xlow   = x - fgConvolutionNSigma * sigma1;
+  Double_t xhigh  = x - fgConvolutionNSigma * sigma1;
+  Double_t step   = (xhigh - xlow) / fgConvolutionSteps;
+  Double_t sum    = 0;
+  
+  for (Int_t i = 0; i <= fgConvolutionSteps/2; i++) { 
+    Double_t x1 = xlow + (i - .5) * step;
+    Double_t x2 = xlow - (i - .5) * step;
+    
+    sum += TMath::Landau(x1, deltap, xi, kTRUE) * TMath::Gaus(x, x1, sigma1);
+    sum += TMath::Landau(x2, deltap, xi, kTRUE) * TMath::Gaus(x, x2, sigma1);
+  }
+  return step * sum * invSq2pi / sigma1;
+}
+
+//____________________________________________________________________
+Double_t 
+AliForwardUtil::NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
+                           Double_t sigma, Double_t sigma_n, Int_t n, 
+                           Double_t* a)
+{
+  Double_t result = LandauGaus(x, delta, xi, sigma, sigma_n);
+  for (Int_t i = 2; i <= n; i++) { 
+    Double_t delta_i =  i * (delta + xi * TMath::Log(i));
+    Double_t xi_i    =  i * xi;
+    Double_t sigma_i =  TMath::Sqrt(Double_t(n)*sigma);
+    Double_t a_i     =  a[i];
+    result           += a_i * AliForwardUtil::LandauGaus(x, delta_i, xi_i, 
+                                                        sigma_i, sigma_n);
+  }
+  return result;
+}
+
+//====================================================================
+AliForwardUtil::ELossFitter::ELossFitter(Double_t lowCut, 
+                                        Double_t maxRange, 
+                                        UShort_t minusBins) 
+  : fLowCut(lowCut), fMaxRange(maxRange), fMinusBins(minusBins), 
+    fFitResults(0), fFunctions(0)
+{
+  fFitResults.SetOwner();
+  fFunctions.SetOwner();
+}
+//____________________________________________________________________
+AliForwardUtil::ELossFitter::~ELossFitter()
+{
+  fFitResults.Delete();
+  fFunctions.Delete();
+}
+//____________________________________________________________________
+void
+AliForwardUtil::ELossFitter::Clear()
+{
+  fFitResults.Clear();
+  fFunctions.Clear();
+}
+//____________________________________________________________________
+TF1*
+AliForwardUtil::ELossFitter::Fit1Particle(TH1* dist, Double_t sigman)
+{
+  // Clear the cache 
+  Clear();
+  
+  // Find the fit range 
+  dist->GetXaxis()->SetRangeUser(fLowCut, fMaxRange);
+  
+  // Normalize peak to 1 
+  Double_t max = dist->GetMaximum(); 
+  dist->Scale(1/max);
+  
+  // Get the bin with maximum 
+  Int_t    maxBin = dist->GetMaximumBin();
+  Double_t maxE   = dist->GetBinLowEdge(maxBin);
+  
+  // Get the low edge 
+  dist->GetXaxis()->SetRangeUser(fLowCut, maxE);
+  Int_t    minBin = maxBin - fMinusBins; // dist->GetMinimumBin();
+  Double_t minE   = TMath::Max(dist->GetBinCenter(minBin),fLowCut);
+  Double_t maxEE  = dist->GetBinCenter(maxBin+2*fMinusBins);
+
+  // Restore the range 
+  dist->GetXaxis()->SetRangeUser(0, fMaxRange);
+  
+  // Define the function to fit 
+  TF1*          landau1 = new TF1("landau1", landauGaus1, minE, maxEE, 5);
+
+  // Set initial guesses, parameter names, and limits  
+  landau1->SetParameters(5,.5,.07,1,sigman);
+  landau1->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}");
+  landau1->SetParLimits(1, minE, fMaxRange);
+  landau1->SetParLimits(2, 0.00, fMaxRange);
+  landau1->SetParLimits(3, 0.01, fMaxRange);
+  if (sigman <= 0)  landau1->FixParameter(4, 0);
+  else              landau1->SetParLimits(4, 0, fMaxRange);
+
+  // Do the fit, getting the result object 
+  TFitResultPtr r = dist->Fit(landau1, "RNQS", "", minE, maxEE);
+
+  fFitResults.AddAtAndExpand(new TFitResult(*r), 0);
+  fFunctions.AddAtAndExpand(landau1, 0);
+
+  return landau1;
+}
+//____________________________________________________________________
+TF1*
+AliForwardUtil::ELossFitter::FitNParticle(TH1* dist, UShort_t n, 
+                                         Double_t sigman)
+{
+  // Get the seed fit result 
+  TFitResult* r = static_cast<TFitResult*>(fFitResults.At(0));
+  TF1*        f = static_cast<TF1*>(fFunctions.At(0));
+  if (!r || !f) { 
+    f = Fit1Particle(dist, sigman);
+    r = static_cast<TFitResult*>(fFitResults.At(0));
+    if (!r || !f) { 
+      ::Warning("FitNLandau", "No first shot at landau fit");
+      return 0;
+    }
+  }
+
+  // Get some parameters from seed fit 
+  Double_t delta1  = r->Parameter(1);
+  Double_t xi1     = r->Parameter(2);
+  Double_t maxEi   = n * (delta1 + xi1 * TMath::Log(n)) + 2 * n * xi1;
+  Double_t minE    = f->GetXmin();
+
+  // Make the fit function 
+  TF1* landaun     = new TF1(Form("landau%d", n), &landauGausN,minE,maxEi,5+n);
+  landaun->SetLineStyle((n % 10)+1);
+  landaun->SetLineWidth(2);
+  landaun->SetParNames("C","#Delta_{p}","#xi", "#sigma", "#sigma_{n}", "N");
+
+  // Set the initial parameters from the seed fit 
+  landaun->SetParameter(0, r->Parameter(0)); // Constant
+  landaun->SetParameter(1, r->Parameter(1)); // Delta 
+  landaun->SetParameter(2, r->Parameter(2)); // xi
+  landaun->SetParameter(3, r->Parameter(3)); // sigma
+  landaun->SetParameter(4, r->Parameter(4)); // sigma_n
+  landaun->SetParLimits(1, minE, fMaxRange); // Delta
+  landaun->SetParLimits(2, 0.00, fMaxRange); // xi
+  landaun->SetParLimits(3, 0.01, fMaxRange); // sigma
+  landaun->SetParLimits(4, 0.00, fMaxRange); // sigma_n
+  // Fix the number parameter 
+  landaun->FixParameter(5, n);               // N
+
+  // Set the range and name of the scale parameters 
+  for (UShort_t i = 2; i <= n; i++) {// Take parameters from last fit 
+    landaun->SetParameter(5+i-2, n == 2 ? 0.05 : 0);
+    landaun->SetParLimits(5+i-2, 0,1);
+    landaun->SetParName(5+i-2, Form("a_{%d}", i));
+  }
+
+  // Do the fit 
+  TFitResultPtr tr = dist->Fit(landaun, "RSQN", "", minE, maxEi);
+  
+  fFitResults.AddAtAndExpand(new TFitResult(*tr), n-1);
+  fFunctions.AddAtAndExpand(landaun, n-1);
+  
+  return landaun;
+}  
 
 //====================================================================
 AliForwardUtil::Histos::~Histos()
@@ -71,19 +297,6 @@ AliForwardUtil::Histos::Get(UShort_t d, Char_t r) const
   }
   return 0;
 }
-#if 0
-//____________________________________________________________________
-TH2D*
-AliForwardUtil::Histos::Get(UShort_t d, Char_t r)
-{
-  switch (d) { 
-  case 1: return fFMD1i;
-  case 2: return (r == 'I' || r == 'i' ? fFMD2i : fFMD2o);
-  case 3: return (r == 'I' || r == 'i' ? fFMD3i : fFMD3o);
-  }
-  default: return 0;
-}
-#endif
 //====================================================================
 TList*
 AliForwardUtil::RingHistos::DefineOutputList(TList* d) const
@@ -110,138 +323,6 @@ AliForwardUtil::RingHistos::GetOutputHist(TList* d, const char* name) const
   return static_cast<TH1*>(d->FindObject(name));
 }
 
-//====================================================================
-Bool_t
-AliForwardUtil::ReadTriggers(AliESDEvent* esd, UInt_t& triggers,
-                            TH1I* hTriggers)
-{
-  triggers = 0;
-
-  // Get the analysis manager - should always be there 
-  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
-  if (!am) { 
-    AliWarningGeneral("ReadTriggers","No analysis manager defined!");
-    return kFALSE;
-  }
-
-  // Get the input handler - should always be there 
-  AliInputEventHandler* ih = 
-    static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
-  if (!ih) { 
-    AliWarningGeneral("ReadTriggers","No input handler");
-    return kFALSE;
-  }
-  
-  // Get the physics selection - add that by using the macro 
-  // AddTaskPhysicsSelection.C 
-  AliPhysicsSelection* ps = 
-    static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
-  if (!ps) { 
-    AliWarningGeneral("ReadTriggers","No physics selection");
-    return kFALSE;
-  }
-  
-  // Check if this is a collision candidate (INEL)
-  Bool_t inel = ps->IsCollisionCandidate(esd);
-  if (inel) { 
-    triggers |= AliAODForwardMult::kInel;
-    hTriggers->Fill(0.5);
-  }
-
-  // IF this is inel, see if we have a tracklet 
-  if (inel) { 
-    const AliMultiplicity* spdmult = esd->GetMultiplicity();
-    if (!spdmult) {
-      AliWarningGeneral("ReadTriggers","No SPD multiplicity");
-    }
-    else { 
-      Int_t n = spdmult->GetNumberOfTracklets();
-      for (Int_t j = 0; j < n; j++) { 
-       if(TMath::Abs(spdmult->GetEta(j)) < 1) { 
-         triggers |= AliAODForwardMult::kInelGt0;
-         hTriggers->Fill(1.5);
-         break;
-       }
-      }
-    }
-  }
-
-  // Analyse some trigger stuff 
-  AliTriggerAnalysis ta;
-  if (ta.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kNSD1)) {
-    triggers |= AliAODForwardMult::kNSD;
-    hTriggers->Fill(2.5);
-  }
-
-  // Get trigger stuff 
-  TString trigStr = esd->GetFiredTriggerClasses();
-  if (trigStr.Contains("CBEAMB-ABCE-NOPF-ALL")) {
-    triggers |= AliAODForwardMult::kEmpty;
-    hTriggers->Fill(3.5);
-  }
-
-  if (trigStr.Contains("CINT1A-ABCE-NOPF-ALL")) {
-    triggers |= AliAODForwardMult::kA;
-    hTriggers->Fill(4.5);
-  }
-
-  if (trigStr.Contains("CINT1B-ABCE-NOPF-ALL")) {
-    triggers |= AliAODForwardMult::kB;
-    hTriggers->Fill(5.5);
-  }
-
-
-  if (trigStr.Contains("CINT1C-ABCE-NOPF-ALL")) {
-    triggers |= AliAODForwardMult::kC;
-    hTriggers->Fill(6.5);
-  }
-
-  if (trigStr.Contains("CINT1-E-NOPF-ALL")) {
-    triggers |= AliAODForwardMult::kE;
-    hTriggers->Fill(7.5);
-  }
-
-  return kTRUE;
-}
-//____________________________________________________________________
-Bool_t
-AliForwardUtil::ReadVertex(AliESDEvent* esd, Double_t& vz, Double_t maxErr)
-{
-  // Get the vertex 
-  const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
-  if (!vertex) { 
-#ifdef VERBOSE
-    AliWarningGeneral("ReadVertex","No SPD vertex found in ESD");
-#endif
-    return kFALSE;
-  }
-
-  // Check that enough tracklets contributed 
-  if(vertex->GetNContributors() <= 0) {
-#ifdef VERBOSE
-    AliWarningGeneral("ReadVertex",
-                     Form("Number of contributors to vertex is %d<=0",
-                          vertex->GetNContributors()));
-#endif
-    return kFALSE;
-  }
-
-  // Check that the uncertainty isn't too large 
-  if (vertex->GetZRes() > maxErr) { 
-#ifdef VERBOSE
-    AliWarningGeneral("ReadVertex",
-                     Form("Uncertaintity in Z of vertex is too large %f > %d", 
-                          vertex->GetZRes(), maxErr));
-#endif
-    return kFALSE;
-  }
-
-  // Get the z coordiante 
-  vz = vertex->GetZ();
-              
-  return kTRUE;
-}
-
 //
 // EOF
 //
index 996704e1d84f46a18736fc2bc6f540282381bdb1..d0a7c39cd48ae4c0e4d095e9bf7ce672af4d4fc7 100644 (file)
@@ -2,9 +2,11 @@
 #define ALIROOT_PWG2_FORWARD_ALIFORWARDUTIL_H
 #include <TObject.h>
 #include <TString.h>
+#include <TObjArray.h>
 class TH2D;
 class TH1I;
 class TH1;
+class TF1;
 class TAxis;
 class AliESDEvent;
 
@@ -16,27 +18,160 @@ class AliESDEvent;
 class AliForwardUtil : public TObject
 {
 public:
+  //__________________________________________________________________
+  /**
+   * Number of steps to do in the Landau, Gaussiam convolution 
+   */
+  static Int_t fgConvolutionSteps;
+  //------------------------------------------------------------------
+  /** 
+   * How many sigma's of the Gaussian in the Landau, Gaussian
+   * convolution to integrate over
+   */
+  static Double_t fgConvolutionNSigma;
+  //------------------------------------------------------------------
+  /** 
+   * Calculate the shifted Landau
+   * @f[
+   *    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
+   * @f]
+   *
+   * where @f$ f_{L}@f$ is the ROOT implementation of the Landau
+   * distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
+   * @f$\Delta=0,\xi=1@f$. 
+   *
+   * @param x      Where to evaluate @f$ f'_{L}@f$ 
+   * @param delta  Most probable value 
+   * @param xi     The 'width' of the distribution 
+   *
+   * @return @f$ f'_{L}(x;\Delta,\xi) 
+   */
+  static Double_t Landau(Double_t x, Double_t delta, Double_t xi);
+  
+  //------------------------------------------------------------------
   /** 
-   * Read the trigger information from the ESD event 
+   * Calculate the value of a Landau convolved with a Gaussian 
    * 
-   * @param esd        ESD event 
-   * @param triggers   On return, contains the trigger bits 
-   * @param hTriggers  Histogram to fill 
+   * @f[ 
+   *  f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
+   *    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
+   *                       \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
+   * @f]
    * 
-   * @return @c true on success, @c false otherwise 
+   * where @f$ f'_{L}@ is the Landau distribution, @f$ \Delta@f$ the
+   * energy loss, @f$ \xi@f the width of the Landau, and @f$
+   * \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
+   * variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling 
+   * noise in the detector.  
+   *
+   * Note that this function uses the constants fgConvolutionSteps and
+   * fgConvolutionNSigma
+   * 
+   * References: 
+   *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
+   *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
+   *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
+   * 
+   * @param x         where to evaluate @f$ f@f$
+   * @param delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
+   * @param xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
+   * @param sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f
+   * @param sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f
+   * 
+   * @return @f$ f@f$ evaluated at @f$ x@f$.  
    */
-  static Bool_t ReadTriggers(AliESDEvent* esd, UInt_t& triggers, 
-                            TH1I* hTriggers);
+  static Double_t LandauGaus(Double_t x, Double_t delta, Double_t xi, 
+                            Double_t sigma, Double_t sigma_n);
+  
+  //------------------------------------------------------------------
   /** 
-   * Read the vertex information from the ESD event 
+   * Evaluate 
+   * @f$ 
+   *     f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f(x;\Delta_i,\xi_i,\sigma'_i)
+   * @f$ 
    * 
-   * @param esd  ESD event 
-   * @param vz   On return, the vertex Z position 
+   * where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
+   * Landau with a Gaussian (see LandauGaus).  Note that 
+   * @f$ a_1 = 1@f$, @\Delta_i = i(\Delta_1 + \xi\log(i))@f$, 
+   * @f$\xi_i=i\xi_1@f, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$. 
+   *  
+   * References: 
+   *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
+   *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
+   *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
    * 
-   * @return @c true on success, @c false otherwise 
+   * @param x        Where to evaluate @f$ f_N@f$
+   * @param delta    @f$ \Delta_1@f$ 
+   * @param xi       @f$ \xi_1@f$
+   * @param sigma    @f$ \sigma_1@f$ 
+   * @param sigma_n  @f$ \sigma_n@f$ 
+   * @param n        @f$ N@f$ in the sum above.
+   * @param a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
+   *                 @f$ i > 1@f$ 
+   * 
+   * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$ 
    */
-  static Bool_t ReadVertex(AliESDEvent* esd, Double_t& vz,
-                          Double_t maxErr=0.1);
+  static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
+                             Double_t sigma, Double_t sigma_n, Int_t n, 
+                             Double_t* a);
+  //__________________________________________________________________
+  /** 
+   * Structure to do fits to the energy loss spectrum 
+   * 
+   */
+  struct ELossFitter 
+  {
+    /** 
+     * Constructor 
+     * 
+     * @param lowCut     Lower cut of spectrum - data below this cuts is ignored
+     * @param maxRange   Maximum range to fit to 
+     * @param minusBins  The number of bins below maximum to use 
+     */
+    ELossFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins); 
+    virtual ~ELossFitter();
+    /** 
+     * Clear internal arrays 
+     * 
+     */
+    void Clear();
+    /** 
+     * Fit a 1-particle signal to the passed energy loss distribution 
+     * 
+     * Note that this function clears the internal arrays first 
+     *
+     * @param dist    Data to fit the function to 
+     * @param sigman If larger than zero, the initial guess of the
+     *               detector induced noise. If zero or less, then this 
+     *               parameter is ignored in the fit (fixed at 0)
+     * 
+     * @return The function fitted to the data 
+     */
+    TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
+    /** 
+     * Fit a N-particle signal to the passed energy loss distribution 
+     *
+     * If there's no 1-particle fit present, it does that first 
+     *
+     * @param dist   Data to fit the function to 
+     * @param N      Number of particle signals to fit 
+     * @param sigman If larger than zero, the initial guess of the
+     *               detector induced noise. If zero or less, then this 
+     *               parameter is ignored in the fit (fixed at 0)
+     * 
+     * @return The function fitted to the data 
+     */
+    TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
+     
+
+    const Double_t fLowCut;     // Lower cut on data 
+    const Double_t fMaxRange;   // Maximum range to fit 
+    const UShort_t fMinusBins;  // Number of bins from maximum to fit 1st peak
+    TObjArray fFitResults;      // Array of fit results 
+    TObjArray fFunctions;       // Array of functions 
+  };
+      
+
   //__________________________________________________________________
   /** 
    * Structure to hold histograms 
@@ -137,6 +272,12 @@ public:
     TList* DefineOutputList(TList* d) const;
     TList* GetOutputList(TList* d) const;
     TH1* GetOutputHist(TList* d, const char* name) const;
+    Color_t Color() const 
+    { 
+      return ((fDet == 1 ? kRed : (fDet == 2 ? kGreen : kBlue))
+             + ((fRing == 'I' || fRing == 'i') ? 2 : -2));
+    }
+
     UShort_t fDet; 
     Char_t   fRing;
     TString  fName;