21-jul-2006 NvE AliSample extended for computation of median.
authornick <nick@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 28 Jul 2006 13:35:54 +0000 (13:35 +0000)
committernick <nick@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 28 Jul 2006 13:35:54 +0000 (13:35 +0000)
07-jul-2006 NvE Track Id correctly set in IceLinefit.
28-jul-2006 NvE New class IceChi2 introduced and example macro icechi2.cc added.
                Also IcePandel updated to use a Convoluted Pandel pdf.

12 files changed:
RALICE/AliSample.cxx
RALICE/AliSample.h
RALICE/history.txt
RALICE/icepack/ICEHeaders.h
RALICE/icepack/ICELinkDef.h
RALICE/icepack/IceChi2.cxx [new file with mode: 0644]
RALICE/icepack/IceChi2.h [new file with mode: 0644]
RALICE/icepack/IceLinefit.cxx
RALICE/icepack/IcePandel.cxx
RALICE/icepack/IcePandel.h
RALICE/icepack/history.txt
RALICE/icepack/macros/icechi2.cc [new file with mode: 0644]

index 62aa9c4..5a1abc8 100644 (file)
@@ -56,18 +56,43 @@ AliSample::AliSample()
  fNames[1]='Y';
  fNames[2]='Z';
  fN=0;
+ fStore=0;
+ fX=0;
+ fY=0;
+ fZ=0;
+ fArr=0;
  Reset();
 }
 ///////////////////////////////////////////////////////////////////////////
 AliSample::~AliSample()
 {
 // Default destructor
+ if (fX)
+ {
+  delete fX;
+  fX=0;
+ }
+ if (fY)
+ {
+  delete fY;
+  fY=0;
+ }
+ if (fZ)
+ {
+  delete fZ;
+  fZ=0;
+ }
+ if (fArr)
+ {
+  delete fArr;
+  fArr=0;
+ }
 }
 ///////////////////////////////////////////////////////////////////////////
 void AliSample::Reset()
 {
 // Resetting the statistics values for a certain Sample object
-// Dimension is NOT changed
+// Dimension and storage flag are NOT changed
  fN=0;
  for (Int_t i=0; i<fDim; i++)
  {
@@ -82,6 +107,11 @@ void AliSample::Reset()
    fCor[i][j]=0.;
   }
  }
+
+ // Set storage arrays to initial size
+ if (fX) fX->Set(10);
+ if (fY) fY->Set(10);
+ if (fZ) fZ->Set(10);
 }
 ///////////////////////////////////////////////////////////////////////////
 void AliSample::Enter(Float_t x)
@@ -104,6 +134,16 @@ void AliSample::Enter(Float_t x)
   fN+=1;
   fSum[0]+=x;
   fSum2[0][0]+=x*x;
+
+  // Store all entered data when storage mode has been selected 
+  if (fStore)
+  {
+   if (!fX) fX=new TArrayF(10);
+   if (fX->GetSize() < fN) fX->Set(fN+10);
+   fX->AddAt(x,fN-1);
+  }
+
+  // Compute the various statistics
   Compute();
  }
 }
@@ -111,6 +151,9 @@ void AliSample::Enter(Float_t x)
 void AliSample::Remove(Float_t x)
 {
 // Removing a value from a 1-dim. sample
+
+ if (!fN) return;
+
  if (fDim != 1)
  {
   cout << " *AliSample::remove* Error : Not a 1-dim sample." << endl;
@@ -120,6 +163,25 @@ void AliSample::Remove(Float_t x)
   fN-=1;
   fSum[0]-=x;
   fSum2[0][0]-=x*x;
+
+  // Remove data entry from the storage
+  if (fStore && fX)
+  {
+   Float_t val=0;
+   for (Int_t i=0; i<=fN; i++)
+   {
+    val=fX->At(i);
+    if (fabs(x-val)>1.e-10) continue;
+
+    for (Int_t j=i+1; j<=fN; j++)
+    {
+     val=fX->At(j);
+     fX->AddAt(val,j-1);
+    }
+   }
+  }
+
+  // Compute the various statistics
   Compute();
  }
 }
@@ -148,6 +210,19 @@ void AliSample::Enter(Float_t x,Float_t y)
   fSum2[0][1]+=x*y;
   fSum2[1][0]+=y*x;
   fSum2[1][1]+=y*y;
+
+  // Store all entered data when storage mode has been selected 
+  if (fStore)
+  {
+   if (!fX) fX=new TArrayF(10);
+   if (fX->GetSize() < fN) fX->Set(fN+10);
+   fX->AddAt(x,fN-1);
+   if (!fY) fY=new TArrayF(10);
+   if (fY->GetSize() < fN) fY->Set(fN+10);
+   fY->AddAt(y,fN-1);
+  }
+
+  // Compute the various statistics
   Compute();
  }
 }
@@ -155,6 +230,9 @@ void AliSample::Enter(Float_t x,Float_t y)
 void AliSample::Remove(Float_t x,Float_t y)
 {
 // Removing a pair (x,y) from a 2-dim. sample
+
+ if (!fN) return;
+
  if (fDim != 2)
  {
   cout << " *AliSample::remove* Error : Not a 2-dim sample." << endl;
@@ -168,6 +246,29 @@ void AliSample::Remove(Float_t x,Float_t y)
   fSum2[0][1]-=x*y;
   fSum2[1][0]-=y*x;
   fSum2[1][1]-=y*y;
+
+  // Remove data entry from the storage
+  if (fStore && fX && fY)
+  {
+   Float_t val=0;
+   for (Int_t i=0; i<=fN; i++)
+   {
+    val=fX->At(i);
+    if (fabs(x-val)>1.e-10) continue;
+    val=fY->At(i);
+    if (fabs(y-val)>1.e-10) continue;
+
+    for (Int_t j=i+1; j<=fN; j++)
+    {
+     val=fX->At(j);
+     fX->AddAt(val,j-1);
+     val=fY->At(j);
+     fY->AddAt(val,j-1);
+    }
+   }
+  }
+
+  // Compute the various statistics
   Compute();
  }
 }
@@ -202,6 +303,22 @@ void AliSample::Enter(Float_t x,Float_t y,Float_t z)
   fSum2[2][0]+=z*x;
   fSum2[2][1]+=z*y;
   fSum2[2][2]+=z*z;
+
+  // Store all entered data when storage mode has been selected 
+  if (fStore)
+  {
+   if (!fX) fX=new TArrayF(10);
+   if (fX->GetSize() < fN) fX->Set(fN+10);
+   fX->AddAt(x,fN-1);
+   if (!fY) fY=new TArrayF(10);
+   if (fY->GetSize() < fN) fY->Set(fN+10);
+   fY->AddAt(y,fN-1);
+   if (!fZ) fZ=new TArrayF(10);
+   if (fZ->GetSize() < fN) fZ->Set(fN+10);
+   fZ->AddAt(z,fN-1);
+  }
+
+  // Compute the various statistics
   Compute();
  }
 }
@@ -209,6 +326,9 @@ void AliSample::Enter(Float_t x,Float_t y,Float_t z)
 void AliSample::Remove(Float_t x,Float_t y,Float_t z)
 {
 // Removing a set (x,y,z) from a 3-dim. sample
+
+ if (!fN) return;
+
  if (fDim != 3)
  {
   cout << " *AliSample::remove* Error : Not a 3-dim sample." << endl;
@@ -228,6 +348,33 @@ void AliSample::Remove(Float_t x,Float_t y,Float_t z)
   fSum2[2][0]-=z*x;
   fSum2[2][1]-=z*y;
   fSum2[2][2]-=z*z;
+
+  // Remove data entry from the storage
+  if (fStore && fX && fY && fZ)
+  {
+   Float_t val=0;
+   for (Int_t i=0; i<=fN; i++)
+   {
+    val=fX->At(i);
+    if (fabs(x-val)>1.e-10) continue;
+    val=fY->At(i);
+    if (fabs(y-val)>1.e-10) continue;
+    val=fZ->At(i);
+    if (fabs(z-val)>1.e-10) continue;
+
+    for (Int_t j=i+1; j<=fN; j++)
+    {
+     val=fX->At(j);
+     fX->AddAt(val,j-1);
+     val=fY->At(j);
+     fY->AddAt(val,j-1);
+     val=fZ->At(j);
+     fZ->AddAt(val,j-1);
+    }
+   }
+  }
+
+  // Compute the various statistics
   Compute();
  }
 }
@@ -355,18 +502,20 @@ Float_t AliSample::GetCor(Int_t i,Int_t j) const
  }
 }
 ///////////////////////////////////////////////////////////////////////////
-void AliSample::Data() const
+void AliSample::Data()
 {
 // Printing of statistics of all variables
  for (Int_t i=0; i<fDim; i++)
  {
  cout << " " << fNames[i] << " : N = " << fN;
  cout << " Sum = " << fSum[i] << " Mean = " << fMean[i];
- cout << " Var = " << fVar[i] << " Sigma = " << fSigma[i] << endl;
+ cout << " Var = " << fVar[i] << " Sigma = " << fSigma[i];
+ if (fStore) cout << " Median = " << GetMedian(i+1);
+ cout << endl;
  }
 }
 ///////////////////////////////////////////////////////////////////////////
-void AliSample::Data(Int_t i) const
+void AliSample::Data(Int_t i)
 {
 // Printing of statistics of ith variable
  if (fDim < i)
@@ -377,7 +526,9 @@ void AliSample::Data(Int_t i) const
  {
   cout << " " << fNames[i-1] << " : N = " << fN;
   cout << " Sum = " << fSum[i-1] << " Mean = " << fMean[i-1];
-  cout << " Var = " << fVar[i-1] << " Sigma = " << fSigma[i-1] << endl;
+  cout << " Var = " << fVar[i-1] << " Sigma = " << fSigma[i-1];
+  if (fStore) cout << " Median = " << GetMedian(i);
+  cout << endl;
  }
 }
 ///////////////////////////////////////////////////////////////////////////
@@ -397,3 +548,114 @@ void AliSample::Data(Int_t i,Int_t j) const
  }
 }
 ///////////////////////////////////////////////////////////////////////////
+void AliSample::SetStoreMode(Int_t mode)
+{
+// Set storage mode for all entered data.
+//
+// mode = 0 : Entered data will not be stored
+//        1 : All data will be stored as entered
+//
+// By default the storage mode is set to 0 in the constructor of this class.
+// The default at invokation of this memberfunction is mode=1.
+//
+// For normal statistics evaluation (e.g. mean, sigma, covariance etc...)
+// storage of entered data is not needed. This is the default mode
+// of operation and is the most efficient w.r.t. cpu time and memory.
+// However, when calculation of a median is required, then the data
+// storage mode has be activated.
+//
+// Note : Activation of storage mode can only be performed before the
+//        first data item is entered. 
+
+ if (fN)
+ {
+  cout << " *AliSample::SetStore* Storage mode can only be set before first data." << endl;
+ }
+ else
+ {
+  if (mode==0 || mode==1) fStore=mode;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Int_t AliSample::GetStoreMode() const
+{
+// Provide the storage mode
+ return fStore;
+}
+///////////////////////////////////////////////////////////////////////////
+Float_t AliSample::GetMedian(Int_t i)
+{
+// Provide the median of a certain variable
+
+ if (fDim < i)
+ {
+  cout << " *AliSample::mean* Error : Dimension less than " << i << endl;
+  return 0;
+ }
+
+ if (!fStore)
+ {
+  cout << " *AliSample::GetMedian* Error : Storage of data entries was not activated." << endl;
+  return 0;
+ }
+
+ // Prepare temp. array to hold the ordered values
+ if (!fArr)
+ {
+  fArr=new TArrayF(fN);
+ }
+ else
+ {
+  if (fArr->GetSize() < fN) fArr->Set(fN);
+ }
+
+ // Order the values of the specified variable
+ Float_t val=0;
+ Int_t iadd=0;
+ for (Int_t j=0; j<fN; j++)
+ {
+  if (i==1) val=fX->At(j);
+  if (i==2) val=fY->At(j);
+  if (i==3) val=fZ->At(j);
+
+  iadd=0;
+  if (j==0)
+  {
+   fArr->AddAt(val,j);
+   iadd=1;
+  }
+  else
+  {
+   for (Int_t k=0; k<j; k++)
+   {
+    if (val>=fArr->At(k)) continue;
+    // Put value in between the existing ones
+    for (Int_t m=j-1; m>=k; m--)
+    {
+     fArr->AddAt(fArr->At(m),m+1);
+    }
+    fArr->AddAt(val,k);
+    iadd=1;
+    break;
+   }
+
+   if (!iadd)
+   {
+    fArr->AddAt(val,j);
+   }
+  }
+ }
+
+ Float_t median=0;
+ Int_t index=fN/2;
+ if (fN%2) // Odd number of entries
+ {
+  median=fArr->At(index);
+ }
+ else // Even number of entries
+ {
+  median=(fArr->At(index-1)+fArr->At(index))/2.;
+ }
+ return median;
+}
+///////////////////////////////////////////////////////////////////////////
index d9c960f..9f1ca69 100644 (file)
@@ -8,6 +8,7 @@
 #include <math.h>
 
 #include "Rtypes.h"
+#include "TArrayF.h"
  
 class AliSample
 {
@@ -29,9 +30,12 @@ class AliSample
   Float_t GetSigma(Int_t i) const;              // Standard deviation for i-th variable
   Float_t GetCov(Int_t i, Int_t j) const;       // Covariance for i-th and j-th variable
   Float_t GetCor(Int_t i, Int_t j) const;       // Correlation for i-th and j-th variable
-  void Data() const;                            // Stat. info for the complete sample
-  void Data(Int_t i) const;                     // Stat. info for the i-th variable
+  Float_t GetMedian(Int_t i);                   // Provide median for i-th variable
+  void Data();                                  // Stat. info for the complete sample
+  void Data(Int_t i);                           // Stat. info for the i-th variable
   void Data(Int_t i, Int_t j) const;            // Stat. info for i-th and j-th variable
+  void SetStoreMode(Int_t mode=1);              // Set mode for storage of all entered data
+  Int_t GetStoreMode() const;                   // Provide storage mode of all entered data
  
  private:
   Int_t fDim;                      // Dimension of the sample
@@ -46,6 +50,11 @@ class AliSample
   Float_t fSigma[fMaxdim];         // Standard deviation for each variable
   Float_t fCov[fMaxdim][fMaxdim];  // Covariances of pairs of variables
   Float_t fCor[fMaxdim][fMaxdim];  // Correlations of pairs of variables
+  Int_t fStore;                    // Flag to denote storage of all entered data 
+  TArrayF* fX;                     // Storage array for the 1st variable (e.g. X)
+  TArrayF* fY;                     // Storage array for the 2nd variable (e.g. Y)
+  TArrayF* fZ;                     // Storage array for the 3rd variable (e.g. Z)
+  TArrayF* fArr;                   // Temp. storage array for ordering
 
  ClassDef(AliSample,0) // Statistics tools for various multi-dimensional data samples.
 };
index 23a5881..7e229bd 100644 (file)
                 Double_t replaced by Double32_t in Ali3Vector.h, Ali4Vector.h and AliBoost.h
                 and fDresult suppressed for I/O streaming. This will reduce the produced
                 output filesize.
+21-jul-2006 NvE AliSample extended for computation of median.
index 2d95d1e..ac5a431 100644 (file)
@@ -20,4 +20,5 @@
 #include "IceDwalk.h"
 #include "IcePandel.h"
 #include "IceLinefit.h"
+#include "IceChi2.h"
 
index 35d99d4..79f333a 100644 (file)
@@ -25,5 +25,6 @@
  #pragma link C++ class IceDwalk+;
  #pragma link C++ class IcePandel+;
  #pragma link C++ class IceLinefit+;
+ #pragma link C++ class IceChi2+;
 #endif
  
diff --git a/RALICE/icepack/IceChi2.cxx b/RALICE/icepack/IceChi2.cxx
new file mode 100644 (file)
index 0000000..b62ee63
--- /dev/null
@@ -0,0 +1,765 @@
+/*******************************************************************************
+ * Copyright(c) 2003, IceCube Experiment at the South Pole. All rights reserved.
+ *
+ * Author: The IceCube RALICE-based Offline Project.
+ * Contributors are mentioned in the code where appropriate.
+ *
+ * Permission to use, copy, modify and distribute this software and its
+ * documentation strictly for non-commercial purposes is hereby granted
+ * without fee, provided that the above copyright notice appears in all
+ * copies and that both the copyright notice and this permission notice
+ * appear in the supporting documentation.
+ * The authors make no claims about the suitability of this software for
+ * any purpose. It is provided "as is" without express or implied warranty.
+ *******************************************************************************/
+
+// $Id$
+
+///////////////////////////////////////////////////////////////////////////
+// Class IceChi2
+// TTask derived class to perform track fitting via chi-squared minimisation.
+//
+// For the minimisation process the TFitter facility, which is basically Minuit,
+// is used. Minimisation is performed by invokation of the SIMPLEX method,
+// followed by an invokation of HESSE to determine the uncertainties on the results.
+// The statistics of the TFitter result are stored as an AliSignal object
+// in the track, which can be obtained via the GetFitDetails memberfunction.
+// After the chi-squared minimisation procedure has been performed, an overall
+// plausibility for the fitted track will be determined based on a convoluted
+// Pandel pdf value for each used hit.
+// This track plausibility is expressed in terms of a Bayesian psi value
+// w.r.t. a Convoluted Pandel PDF.
+// The Baysian psi value is defined as -loglikelihood in a decibel scale.
+// This implies psi=-10*log10(L) where L=p(D|HI) being the likelihood of
+// the data D under the hypothesis H and prior information I.
+// Since all (associated) hits contribute independently to the Bayesian psi
+// value, this psi value is built up by summation of the various hit contributions.
+// As such, the FitDetails entries contain the statistics of all the different
+// hit contributions, like PsiMedian, PsiMean, and PsiSigma.
+// The Bayesian psi value is available in the fit details under the name "PsiSum".
+// In addition the standard Minuit results like IERFIT, FCN, EDM etc... are
+// also available from the FitDetails.
+//
+// The convoluted Pandel value is evaluated in various areas in the distance-time
+// space as described in the CPandel writeup by O. Fadiran, G. Japaridze
+// and N. van Eijndhoven.
+// In case the distance-time point of a certain hit falls outside the
+// validity rectangle, the point is moved onto the corresponding side location
+// of the rectangle. For this new location the Pandel value is evaluated for
+// this hit and an extra penalty is added to the corresponding psi value
+// for this hit. 
+// By default this penalty value amounts to 0 dB, but the user can
+// modify this penalty value via the memberfunction SetPenalty.
+// This allows investigation/tuning of the sensitivity to hits with
+// extreme distance and/or time residual values.
+//
+// Use the UseTracks memberfunction to specify the first guess tracks
+// to be processed by the minimiser.
+// By default only the first encountered IceDwalk track will be processed.
+//
+// Use the SelectHits memberfunction to specify the hits to be used.
+// By default only the hits associated to the first guess track are used.
+//
+// The fit processor printlevel can be selected via the memberfunction
+// SetPrintLevel. By default all printout is suppressed (i.e. level=-2).
+//
+// An example of how to invoke this processor after Xtalk, hit cleaning
+// and a direct walk first guess estimate can be found in the ROOT example
+// macro icechi2.cc which resides in the /macros subdirectory.
+// 
+// The minimisation results are stored in the IceEvent structure as
+// tracks with as default the name "IceChi2" (just like the first guess
+// results of e.g. IceDwalk).
+// This track name identifier can be modified by the user via the
+// SetTrackName() memberfunction. This will allow unique identification
+// of tracks which are produced when re-processing existing data with
+// different criteria.
+// By default the charge of the produced tracks is set to 0, since
+// no distinction can be made between positive or negative tracks.
+// However, the user can define the track charge by invokation
+// of the memberfunction SetCharge().
+// This facility may be used to distinguish tracks produced by the
+// various reconstruction algorithms in a (3D) colour display
+// (see the class AliHelix for further details).  
+// A pointer to the first guess track which was used as input is available
+// via the GetParentTrack facility of these "IceChi2" tracks.
+// Furthermore, all the hits that were used in the minisation are available
+// via the GetSignal facility of a certain track.
+//
+// An example of how the various data can be accessed is given below,
+// where "evt" indicates the pointer to the IceEvent structure.
+//
+// Example for accessing data :
+// ----------------------------
+// TObjArray* tracks=evt->GetTracks("IceChi2");
+// if (!tracks) return;
+// AliPosition* r0=0;
+// Float_t psi=0;
+// for (Int_t jtk=0; jtk<tracks->GetEntries(); jtk++)
+// {
+//  AliTrack* tx=(AliTrack*)tracks->At(jtk);
+//  if (!tx) continue;
+//  tx->Data("sph");
+//  r0=tx->GetReferencePoint();
+//  if (r0) r0->Data();
+//  sx=(AliSignal*)tx->GetFitDetails();
+//  if (sx) psi=sx->GetSignal("PsiSum");
+//  AliTrack* tx2=tx->GetParentTrack();
+//  if (!tx2) continue;
+//  tx2->Data("sph");
+//  r0=tx2->GetReferencePoint();
+//  if (r0) r0->Data();
+// }
+//
+// Notes :
+// -------
+// 1) This processor only works properly on data which are Time and ADC
+//    calibrated and contain tracks from first guess algorithms like
+//    e.g. IceDwalk.
+// 2) In view of the usage of TFitter/Minuit minimisation, a global pointer
+//    to the instance of this class (gIceChi2) and a global static
+//    wrapper function (IceChi2FCN) have been introduced, to allow the
+//    actual minimisation to be performed via the memberfunction FitFCN.
+//    This implies that in a certain processing job only 1 instance of
+//    this IceChi2 class may occur.
+//
+//--- Author: Nick van Eijndhoven 16-may-2006 Utrecht University
+//- Modified: NvE $Date$ Utrecht University
+///////////////////////////////////////////////////////////////////////////
+#include "IceChi2.h"
+#include "Riostream.h"
+
+// Global pointer to the instance of this object
+ IceChi2* gIceChi2=0;
+
+// TFitter/Minuit interface to IceChi2::FitFCN
+ void IceChi2FCN(Int_t& npar,Double_t* gin,Double_t& f,Double_t* u,Int_t flag)
+ {
+  if (gIceChi2) gIceChi2->FitFCN(npar,gin,f,u,flag);
+ }
+
+ClassImp(IceChi2) // Class implementation to enable ROOT I/O
+
+IceChi2::IceChi2(const char* name,const char* title) : TTask(name,title)
+{
+// Default constructor.
+ fFirst=1;
+ fPrint=-2;
+ fSelhits=1;
+ fEvt=0;
+ fUseNames=0;
+ fUseNtk=0;
+ fHits=0;
+ fFitter=0;
+ fTrackname="IceChi2";
+ fCharge=0;
+ fPenalty=0;
+ fTkfit=0;
+ fFitstats=0;
+
+ // Set the global pointer to this instance
+ gIceChi2=this;
+}
+///////////////////////////////////////////////////////////////////////////
+IceChi2::~IceChi2()
+{
+// Default destructor.
+ if (fUseNames)
+ {
+  delete fUseNames;
+  fUseNames=0;
+ }
+ if (fUseNtk)
+ {
+  delete fUseNtk;
+  fUseNtk=0;
+ }
+ if (fHits)
+ {
+  delete fHits;
+  fHits=0;
+ }
+ if (fFitter)
+ {
+  delete fFitter;
+  fFitter=0;
+ }
+ if (fTkfit)
+ {
+  delete fTkfit;
+  fTkfit=0;
+ }
+ if (fFitstats)
+ {
+  delete fFitstats;
+  fFitstats=0;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+void IceChi2::Exec(Option_t* opt)
+{
+// Implementation of the hit fitting procedure.
+
+ TString name=opt;
+ AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
+
+ if (!parent) return;
+
+ fEvt=(IceEvent*)parent->GetObject("IceEvent");
+ if (!fEvt) return;
+
+ if (!fUseNames) UseTracks("IceDwalk",1);
+
+ Int_t nclasses=fUseNames->GetEntries(); // Number of first guess classes to be processed
+ Int_t ntkmax=0; // Max. number of tracks for a certain class
+ TObjString* strx=0;
+ TString str;
+
+ if (fFirst)
+ {
+  cout << " *IceChi2* First guess selections to be processed (-1=all)." << endl;
+  for (Int_t i=0; i<nclasses; i++)
+  {
+   strx=(TObjString*)fUseNames->At(i);
+   if (!strx) continue;
+   str=strx->GetString();
+   ntkmax=fUseNtk->At(i);
+   cout << " Maximally " << ntkmax << " track(s) per event for procedure : " << str.Data() << endl;
+  }
+  cout << " *IceChi2* Hit selection mode : " << fSelhits << endl;
+  cout << " *IceChi2* Penalty value for psi evaluation outside range : " << fPenalty << endl;
+  cout << endl;
+
+  fPsistats.SetStoreMode(1);
+
+  fFirst=0;
+ }
+
+ const Double_t pi=acos(-1.);
+
+ // Initialisation of the minimisation processor
+ Double_t arglist[100];
+ if (!fFitter) fFitter=new TFitter();
+
+ // The number of reconstructed tracks already present in the event
+ Int_t ntkreco=fEvt->GetNtracks(1);
+
+ if (!fHits)
+ {
+  fHits=new TObjArray();
+ }
+ else
+ {
+  fHits->Clear();
+ }
+
+ // If selected, use all the good quality hits of the complete event
+ if (fSelhits==0)
+ {
+  TObjArray* hits=fEvt->GetHits("IceGOM");
+  for (Int_t ih=0; ih<hits->GetEntries(); ih++)
+  {
+   AliSignal* sx=(AliSignal*)hits->At(ih);
+   if (!sx) continue;
+   if (sx->GetDeadValue("ADC") || sx->GetDeadValue("LE") || sx->GetDeadValue("TOT")) continue;
+   fHits->Add(sx);
+  }
+ }
+
+ // Track by track processing of the selected first guess classes
+ Float_t vec[3];
+ Float_t err[3];
+ Ali3Vector p;
+ AliPosition pos;
+ if (!fTkfit)
+ {
+  fTkfit=new AliTrack();
+  fTkfit->SetNameTitle(fTrackname.Data(),"IceChi2 fit result");
+ }
+ if (!fFitstats)
+ {
+  fFitstats=new AliSignal();
+  fFitstats->SetNameTitle("Fitstats","TFitter stats for Chi2 fit");
+  fFitstats->SetSlotName("IERFIT",1);
+  fFitstats->SetSlotName("FCN",2);
+  fFitstats->SetSlotName("EDM",3);
+  fFitstats->SetSlotName("NVARS",4);
+  fFitstats->SetSlotName("IERERR",5);
+  fFitstats->SetSlotName("PsiSum",6);
+  fFitstats->SetSlotName("PsiMedian",7);
+  fFitstats->SetSlotName("PsiMean",8);
+  fFitstats->SetSlotName("PsiSigma",9);
+ }
+ Float_t x,y,z,theta,phi,t0;
+ Double_t amin,edm,errdef; // Minimisation stats
+ Int_t ierfit,iererr,nvpar,nparx; // Minimisation stats
+ Double_t psi; // Bayesian psi value for the fitted track w.r.t. a Convoluted Pandel PDF
+ Int_t ntk=0;
+ Int_t nsig=0;
+ AliTrack* track=0;
+ for (Int_t iclass=0; iclass<nclasses; iclass++) // Loop over first guess classes
+ {
+  strx=(TObjString*)fUseNames->At(iclass);
+  if (!strx) continue;
+  str=strx->GetString();
+  ntkmax=fUseNtk->At(iclass);
+  TObjArray* tracks=fEvt->GetTracks(str);
+  ntk=tracks->GetEntries();
+  if (ntkmax>0 && ntk>ntkmax) ntk=ntkmax;
+
+  for (Int_t jtk=0; jtk<ntk; jtk++) // Loop over tracks of a certain class
+  {
+   track=(AliTrack*)tracks->At(jtk);
+   if (!track) continue;
+
+   AliPosition* r0=track->GetReferencePoint();
+   if (!r0) continue;
+
+   AliTimestamp* tt0=r0->GetTimestamp();
+
+   // If selected, use only the first guess track associated hits
+   if (fSelhits==1)
+   {
+    fHits->Clear();
+    nsig=track->GetNsignals();
+    for (Int_t is=1; is<=nsig; is++)
+    {
+     AliSignal* sx=track->GetSignal(is);
+     if (!sx) continue;
+     if (!sx->GetDevice()->InheritsFrom("IceGOM")) continue;
+     if (sx->GetDeadValue("ADC") || sx->GetDeadValue("LE") || sx->GetDeadValue("TOT")) continue;
+     fHits->Add(sx);
+    }
+   }
+   
+   if (!fHits->GetEntries()) continue;
+
+   r0->GetVector(vec,"car");
+   r0->GetErrors(err,"car");
+
+   x=vec[0];
+   y=vec[1];
+   z=vec[2];
+
+   p=track->Get3Momentum();
+   p.GetVector(vec,"sph");
+
+   theta=vec[1];
+   phi=vec[2];
+
+   t0=fEvt->GetDifference(tt0,"ns");
+   
+   // Process this first guess track with its associated hits
+   fFitter->Clear();
+
+   // Set user selected TFitter printout level
+   arglist[0]=fPrint;
+   if (fPrint==-2) arglist[0]=-1;
+   fFitter->ExecuteCommand("SET PRINT",arglist,1);
+   if (fPrint==-2) fFitter->ExecuteCommand("SET NOWARNINGS",arglist,0);
+
+   fFitter->SetFitMethod("chisquare");
+
+   fFitter->SetParameter(0,"r0x",x,0.1,0,0);
+   fFitter->SetParameter(1,"r0y",y,0.1,0,0);
+   fFitter->SetParameter(2,"r0z",z,0.1,0,0);
+   fFitter->SetParameter(3,"theta",theta,0.001,0,pi);
+   fFitter->SetParameter(4,"phi",phi,0.001,0,2.*pi);
+   fFitter->SetParameter(5,"t0",t0,1.,0,32000);
+
+   fFitter->SetFCN(IceChi2FCN);
+
+   fTkfit->Reset();
+
+   arglist[0]=0;
+   ierfit=fFitter->ExecuteCommand("SIMPLEX",arglist,0);
+
+   fFitter->GetStats(amin,edm,errdef,nvpar,nparx);
+   fFitstats->Reset();
+   fFitstats->SetSignal(ierfit,1);
+   fFitstats->SetSignal(amin,2);
+   fFitstats->SetSignal(edm,3);
+   fFitstats->SetSignal(nvpar,4);
+
+   iererr=fFitter->ExecuteCommand("HESSE",arglist,0);
+   fFitstats->SetSignal(iererr,5);
+
+
+   // Resulting parameters after minimisation and error calculation
+   vec[0]=fFitter->GetParameter(0);
+   vec[1]=fFitter->GetParameter(1);
+   vec[2]=fFitter->GetParameter(2);
+   err[0]=fFitter->GetParError(0);
+   err[1]=fFitter->GetParError(1);
+   err[2]=fFitter->GetParError(2);
+   pos.SetPosition(vec,"car");
+   pos.SetPositionErrors(err,"car");
+
+   vec[0]=1.;
+   vec[1]=fFitter->GetParameter(3);
+   vec[2]=fFitter->GetParameter(4);
+   err[0]=0.;
+   err[1]=fFitter->GetParError(3);
+   err[2]=fFitter->GetParError(4);
+   p.SetVector(vec,"sph");
+   p.SetErrors(err,"sph");
+
+   t0=fFitter->GetParameter(5);
+   AliTimestamp t0fit((AliTimestamp)(*fEvt));
+   t0fit.Add(0,0,t0);
+
+   // Enter the fit result as a track in the event structure
+   ntkreco++;
+   fTkfit->SetId(ntkreco);
+   fTkfit->SetCharge(fCharge);
+   fTkfit->SetParentTrack(track);
+   pos.SetTimestamp(t0fit);
+   fTkfit->SetTimestamp(t0fit);
+   fTkfit->SetReferencePoint(pos);
+   fTkfit->Set3Momentum(p);
+   for (Int_t ihit=0; ihit<fHits->GetEntries(); ihit++)
+   {
+    AliSignal* sx=(AliSignal*)fHits->At(ihit);
+    if (sx) fTkfit->AddSignal(*sx);
+   }
+   psi=GetPsi(fTkfit);
+   fFitstats->SetSignal(psi,6);
+   fFitstats->SetSignal(fPsistats.GetMedian(1),7);
+   fFitstats->SetSignal(fPsistats.GetMean(1),8);
+   fFitstats->SetSignal(fPsistats.GetSigma(1),9);
+   fTkfit->SetFitDetails(fFitstats);
+   fEvt->AddTrack(fTkfit);
+  } // End loop over tracks
+ } // End loop over first guess classes
+
+}
+///////////////////////////////////////////////////////////////////////////
+void IceChi2::SetPrintLevel(Int_t level)
+{
+// Set the fitter (Minuit) print level.
+// See the TFitter and TMinuit docs for details.
+//
+// Note : level=-2 suppresses also all fit processor warnings.
+//        
+// The default in the constructor is level=-2. 
+
+ fPrint=level;
+}
+///////////////////////////////////////////////////////////////////////////
+void IceChi2::UseTracks(TString classname,Int_t n)
+{
+// Specification of the first guess tracks to be used.
+//
+// classname : Specifies the first guess algorithm (e.g. "IceDwalk");
+// n : Specifies the max. number of these tracks to be used
+//
+// Note : n<0 will use all the existing tracks of the specified classname
+//
+// The default is n=-1.
+//
+// Consecutive invokations of this memberfunction with different classnames
+// will result in an incremental effect.
+//
+// Example :
+// ---------
+// UseTracks("IceDwalk",5);
+// UseTracks("IceLinefit",2);
+// UseTracks("IceJams");
+//
+// This will use the first 5 IceDwalk, the first 2 IceLinefit and all the
+// IceJams tracks which are encountered in the event structure.
+
+ if (!fUseNames)
+ {
+  fUseNames=new TObjArray();
+  fUseNames->SetOwner();
+ }
+ if (!fUseNtk) fUseNtk=new TArrayI();
+
+ // Check if this classname has already been specified before 
+ TString s;
+ Int_t nen=fUseNames->GetEntries();
+ for (Int_t i=0; i<nen; i++)
+ {
+  TObjString* sx=(TObjString*)fUseNames->At(i);
+  if (!sx) continue;
+  s=sx->GetString();
+  if (s==classname) return;
+ }
+
+ // New classname to be added into the storage
+ if (nen >= fUseNames->GetSize()) fUseNames->Expand(nen+1);
+ if (nen >= fUseNtk->GetSize()) fUseNtk->Set(nen+1);
+ TObjString* name=new TObjString();
+ name->SetString(classname);
+ fUseNames->Add(name);
+ fUseNtk->AddAt(n,nen);
+}
+///////////////////////////////////////////////////////////////////////////
+void IceChi2::SelectHits(Int_t mode)
+{
+// Specification of the hits to be used in the minimisation.
+//
+// mode = 0 : All hit cleaning survived hits of the complete event are used
+//        1 : Only the associated hits are used for each first guess track
+//
+// The default is mode=1.
+
+ if (mode==0 || mode==1) fSelhits=mode;
+}
+///////////////////////////////////////////////////////////////////////////
+void IceChi2::SetTrackName(TString s)
+{
+// Set (alternative) name identifier for the produced tracks.
+// This allows unique identification of (newly) produced pandel tracks
+// in case of re-processing of existing data with different criteria.
+// By default the produced tracks have the name "IceChi2" which is
+// set in the constructor of this class.
+ fTrackname=s;
+}
+///////////////////////////////////////////////////////////////////////////
+void IceChi2::SetCharge(Float_t charge)
+{
+// Set user defined charge for the produced tracks.
+// This allows identification of these tracks on color displays.
+// By default the produced tracks have charge=0 which is set in the
+// constructor of this class.
+ fCharge=charge;
+}
+///////////////////////////////////////////////////////////////////////////
+void IceChi2::SetPenalty(Float_t val)
+{
+// Set user defined psi penalty value (in dB) for distance-time points that
+// fall outside the validity rectangle.
+// This allows investigation/tuning of the sensitivity to hits with
+// extreme distance and/or time residual values.
+// By default the penalty val=0 is set in the constructor of this class.
+ fPenalty=val;
+}
+///////////////////////////////////////////////////////////////////////////
+void IceChi2::FitFCN(Int_t&,Double_t*,Double_t& f,Double_t* x,Int_t)
+{
+// The chi-squared function used for the minimisation process.
+
+ const Float_t c=0.299792;           // Light speed in vacuum in meters per ns
+ const Float_t nice=1.35634;         // Refractive index of ice
+ const Float_t thetac=acos(1./nice); // Cherenkov angle (in radians)
+
+ // Assumed PMT timing jitter in ns
+ const Double_t sigt=10;
+
+ f=0;
+
+ // The new r0 and p vectors and t0 from the minimisation
+ Float_t vec[3];
+ AliPosition r0;
+ vec[0]=x[0];
+ vec[1]=x[1];
+ vec[2]=x[2];
+ r0.SetPosition(vec,"car");
+
+ Ali3Vector p;
+ vec[0]=1;
+ vec[1]=x[3];
+ vec[2]=x[4];
+ p.SetVector(vec,"sph");
+
+ Float_t t0=x[5];
+
+ // Construct a track with the new values from the minimisation
+ fTkfit->SetReferencePoint(r0);
+ fTkfit->Set3Momentum(p);
+
+ Int_t nhits=fHits->GetEntries();
+ AliPosition rhit;
+ Ali3Vector r12;
+ Float_t d,dist,thit,tgeo;
+ Double_t tres,chi2;
+ for (Int_t i=0; i<nhits; i++)
+ {
+  AliSignal* sx=(AliSignal*)fHits->At(i);
+  if (!sx) continue;
+  IceGOM* omx=(IceGOM*)sx->GetDevice();
+  if (!omx) continue;
+  rhit=omx->GetPosition();
+  d=fTkfit->GetDistance(rhit);
+  r12=rhit-r0;
+  dist=p.Dot(r12)+d*tan(thetac);
+  tgeo=t0+dist/c;
+  thit=sx->GetSignal("LE",7);
+  tres=thit-tgeo;
+
+  // Chisquared contribution for this hit
+  chi2=pow(tres/sigt,2);
+
+  f+=chi2;
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+Double_t IceChi2::GetPsi(AliTrack* t)
+{
+// Provide Bayesian psi value for a track w.r.t. a Convoluted Pandel PDF.
+// The Baysian psi value is defined as -loglikelihood in a decibel scale.
+// This implies psi=-10*log10(L) where L=p(D|HI) being the likelihood of
+// the data D under the hypothesis H and prior information I.
+// In case of error or incomplete information a psi value of -1 is returned.
+
+ const Float_t c=0.299792;           // Light speed in vacuum in meters per ns
+ const Float_t nice=1.35634;         // Refractive index of ice
+ const Float_t thetac=acos(1./nice); // Cherenkov angle (in radians)
+ const Float_t lambda=33.3;          // Light scattering length in ice
+ const Float_t labs=98;              // Light absorbtion length in ice
+ const Float_t cice=c/nice;          // Light speed in ice in meters per ns
+ const Float_t tau=557;
+ const Double_t rho=((1./tau)+(cice/labs));
+ const Double_t pi=acos(-1.);
+
+ // Assumed PMT timing jitter in ns
+ const Double_t sigma=10;
+
+ Double_t psi=-1;
+
+ if (!t) return psi;
+
+ // The r0 and p vectors from the track
+ AliPosition* ref=t->GetReferencePoint();
+ Ali3Vector p=t->Get3Momentum();
+
+ if (!ref || p.GetNorm()<=0) return psi;
+
+ // The number of associated hits and t0 of the track
+ Int_t nhits=t->GetNsignals();
+ AliTimestamp* tstamp=ref->GetTimestamp();
+
+ if (!nhits || !tstamp) return psi;
+
+ AliPosition r0=ref->GetPosition();
+ Float_t t0=fEvt->GetDifference(tstamp,"ns");
+
+ AliPosition rhit;
+ Ali3Vector r12;
+ Float_t d,dist,thit,tgeo;
+ Double_t tres,ksi,eta,pandel;
+ Double_t cpandel1,cpandel2,cpandel3,cpandel;
+ Double_t z,k,alpha,beta,phi,n1,n2,n3,u; // Function parameters for regions 3 and 4
+ psi=0;
+ Int_t ier;
+ Double_t psihit=0;
+ fPsistats.Reset();
+ for (Int_t i=1; i<=nhits; i++)
+ {
+  AliSignal* sx=t->GetSignal(i);
+  if (!sx) continue;
+  IceGOM* omx=(IceGOM*)sx->GetDevice();
+  if (!omx) continue;
+  rhit=omx->GetPosition();
+  d=t->GetDistance(rhit);
+  ksi=d/lambda;
+  r12=rhit-r0;
+  dist=p.Dot(r12)+d*tan(thetac);
+  tgeo=t0+dist/c;
+  thit=sx->GetSignal("LE",7);
+  tres=thit-tgeo;
+
+  // The Convoluted Pandel function evaluation
+  // For definitions of functions and validity regions, see the
+  // CPandel writeup of O. Fadiran, G. Japaridze and N. van Eijndhoven
+
+  // Move points which are outside the validity rectangle in the (tres,ksi) space
+  // to the edge of the validity rectangle and signal the use of the penalty
+  ier=0;
+  if (tres<-25.*sigma)
+  {
+   tres=-25.*sigma;
+   ier=1;
+  }
+  if (tres>3500)
+  {
+   tres=3500;
+   ier=1;
+  }
+  if (ksi>50)
+  {
+   ksi=50;
+   ier=1;
+  }
+
+  eta=(rho*sigma)-(tres/sigma);
+
+  if (ksi<=0) // The zero distance (ksi=0) axis
+  {
+   cpandel=exp(-tres*tres/(2.*sigma*sigma))/(sigma*sqrt(2.*pi));
+  }
+  else if (ksi<=5 && tres>=-5.*sigma && tres<=30.*sigma) // The exact expression in region 1
+  {
+   cpandel1=pow(rho,ksi)*pow(sigma,ksi-1.)*exp(-tres*tres/(2.*sigma*sigma))/pow(2.,0.5*(1.+ksi));
+   cpandel2=ROOT::Math::conf_hyperg(ksi/2.,0.5,eta*eta/2.)/TMath::Gamma((ksi+1.)/2.);
+   cpandel3=sqrt(2.)*eta*ROOT::Math::conf_hyperg((ksi+1.)/2.,1.5,eta*eta/2.)/TMath::Gamma(ksi/2.);
+
+   cpandel=cpandel1*(cpandel2-cpandel3);
+  }
+  else if (ksi<=1 && tres>30.*sigma && tres<=3500) // Approximation in region 2
+  {
+   pandel=pow(rho,ksi)*pow(tres,(ksi-1.))*exp(-rho*tres)/TMath::Gamma(ksi);
+
+   cpandel=exp(rho*rho*sigma*sigma/2.)*pandel;
+  }
+  else if (ksi<=1 && tres<-5.*sigma && tres>=-25.*sigma) // Approximation in region 5
+  {
+   cpandel=pow(rho*sigma,ksi)*pow(eta,-ksi)*exp(-tres*tres/(2.*sigma*sigma))/(sigma*sqrt(2.*pi));
+  }
+  else if (ksi<=50 && tres>=0 && tres<=3500) // Approximation in region 3
+  {
+   z=-eta/sqrt(4.*ksi-2.);
+   k=0.5*(z*sqrt(1.+z*z)+log(z+sqrt(1.+z*z)));
+   alpha=-tres*tres/(2.*sigma*sigma)+eta*eta/4.-ksi/2.+0.25+k*(2.*ksi-1.);
+   alpha+=-log(1.+z*z)/4.-ksi*log(2.)/2.+(ksi-1.)*log(2.*ksi-1.)/2.+ksi*log(rho)+(ksi-1.)*log(sigma);
+   beta=0.5*(z/sqrt(1.+z*z)-1.);
+   n1=beta*(20.*beta*beta+30.*beta+9.)/12.;
+   n2=pow(beta,2.)*(6160.*pow(beta,4.)+18480.*pow(beta,3.)+19404.*pow(beta,2.)+8028.*beta+945.)/288.;
+   n3=27227200.*pow(beta,6.)+122522400.*pow(beta,5.)+220540320.*pow(beta,4.);
+   n3+=200166120.*pow(beta,3.)+94064328.*pow(beta,2.)+20546550.*beta+1403325.;
+   n3*=pow(beta,3.)/51840.;
+   phi=1.-n1/(2.*ksi-1.)+n2/pow(2.*ksi-1.,2.)-n3/pow(2.*ksi-1.,3.);
+   cpandel=exp(alpha)*phi/TMath::Gamma(ksi);
+  }
+  else if (ksi<=50 && tres<0 && tres>=-25.*sigma) // Approximation in region 4
+  {
+   z=eta/sqrt(4.*ksi-2.);
+   k=0.5*(z*sqrt(1.+z*z)+log(z+sqrt(1.+z*z)));
+   u=exp(ksi/2.-0.25)*pow(2.*ksi-1.,-ksi/2.)*pow(2.,(ksi-1.)/2.);
+   beta=0.5*(z/sqrt(1.+z*z)-1.);
+   n1=beta*(20.*beta*beta+30.*beta+9.)/12.;
+   n2=pow(beta,2.)*(6160.*pow(beta,4.)+18480.*pow(beta,3.)+19404.*pow(beta,2.)+8028.*beta+945.)/288.;
+   n3=27227200.*pow(beta,6.)+122522400.*pow(beta,5.)+220540320.*pow(beta,4.);
+   n3+=200166120.*pow(beta,3.)+94064328.*pow(beta,2.)+20546550.*beta+1403325.;
+   n3*=pow(beta,3.)/51840.;
+   phi=1.+n1/(2.*ksi-1.)+n2/pow(2.*ksi-1.,2.)+n3/pow(2.*ksi-1.,3.);
+   cpandel=pow(rho,ksi)*pow(sigma,ksi-1.)*exp(-pow(tres,2.)/(2.*pow(sigma,2.))+pow(eta,2.)/4.)/sqrt(2.*pi);
+   cpandel*=u*phi*exp(-k*(2.*ksi-1.))*pow(1.+z*z,-0.25);
+  }
+  else // (tres,ksi) outside validity rectangle
+  {
+   ier=1;
+   cout << " *IceChi2::GetPsi* Outside rectangle. We should never get here." << endl;
+  }
+
+  // Use 10*log10 expression to obtain intuitive dB scale
+  // Omit (small) negative values which are possible due to computer accuracy
+  psihit=0;
+  if (cpandel>0) psihit=-10.*log10(cpandel);
+
+  // Penalty in dB for (tres,ksi) points outside the validity rectangle
+  if (ier) psihit+=fPenalty;
+
+  // Update the psi statistics for this hit
+  fPsistats.Enter(float(psihit));
+  psi+=psihit;
+ }
+ return psi;
+}
+///////////////////////////////////////////////////////////////////////////
diff --git a/RALICE/icepack/IceChi2.h b/RALICE/icepack/IceChi2.h
new file mode 100644 (file)
index 0000000..c02766a
--- /dev/null
@@ -0,0 +1,54 @@
+#ifndef IceChi2_h
+#define IceChi2_h
+
+// Copyright(c) 2003, IceCube Experiment at the South Pole, All rights reserved.
+// See cxx source for full Copyright notice.
+
+// $Id$
+
+#include "TROOT.h"
+#include "TTask.h"
+#include "TString.h"
+#include "TObjString.h"
+#include "TFitter.h"
+#include "Math/SpecFuncMathMore.h"
+
+#include "AliJob.h"
+#include "AliSample.h"
+#include "IceEvent.h"
+#include "IceGOM.h"
+
+class IceChi2 : public TTask
+{
+ public :
+  IceChi2(const char* name="",const char* title=""); // Constructor
+  virtual ~IceChi2();                                // Destructor
+  virtual void Exec(Option_t* opt);                  // Perform the fitting procedure
+  void SetPrintLevel(Int_t level);                   // Set the fitter (Minuit) printlevel
+  void UseTracks(TString classname,Int_t n=-1);      // Specify first guess tracks to be used
+  void SelectHits(Int_t mode=1);                     // Specify which hits to be used
+  void SetTrackName(TString s);  // Set (alternative) name for the produced tracks
+  void SetCharge(Float_t charge);// Set user defined charge for the produced tracks
+  void SetPenalty(Float_t val);  // Set psi penalty value (dB) for extreme distance and/or time values
+  Double_t GetPsi(AliTrack* t);  // Provide Bayesian psi value for a track w.r.t. CPandel PDF
+  void FitFCN(Int_t&,Double_t*,Double_t&,Double_t*,Int_t); // The minimisation FCN
+
+ protected :
+  Int_t fFirst;         // Flag to denote first invokation of the processor
+  Int_t fPrint;         // Flag to denote the fitter (Minuit) printlevel
+  Int_t fSelhits;       // Flag to denote which hits to be used
+  IceEvent* fEvt;       // Pointer to the current event structure
+  TObjArray* fUseNames; // The first guess classnames to be used 
+  TArrayI* fUseNtk;     // The max. numbers of the various first guess tracks to be used
+  TObjArray* fHits;     // The various hits to be used in the fitting process 
+  TFitter* fFitter;     // Pointer to the minimisation processor
+  TString fTrackname;   // The name identifier for the produced tracks
+  Float_t fCharge;      // User defined charge of the produced tracks
+  Float_t fPenalty;     // User defined psi penalty value (dB) for extreme distance and/or time values
+  AliTrack* fTkfit;     // Pointer to the produced fitted track
+  AliSignal* fFitstats; // The fit details of the produced fitted track
+  AliSample fPsistats;  // Statistics of the Bayesian psi value for the best fitted track
+
+ ClassDef(IceChi2,1) // TTask derived class to perform Chi-squared track fitting
+};
+#endif
index 6ee1191..dd38447 100644 (file)
@@ -207,6 +207,8 @@ void IceLinefit::Exec(Option_t* opt)
  AliTrack* trk=evt->GetTrack(evt->GetNtracks());
  if (!trk) return;
 
+ trk->SetId(evt->GetNtracks(1)+1);
+
  Ali3Vector p;
  Float_t vec[3];
  v.GetVector(vec,"sph");
index 19f53ce..9f0699c 100644 (file)
 
 ///////////////////////////////////////////////////////////////////////////
 // Class IcePandel
-// TTask derived class to perform Pandel fitting.
+// TTask derived class to perform track fitting via minimisation of a
+// convoluted Pandel pdf.
 //
 // The code in this processor is based on the algorithms as developed
-// by Oladipo Fadiran and George Japaridze (Clark Atlanta University, USA).
+// by Oladipo Fadiran, George Japaridze (Clark Atlanta University, USA)
+// and Nick van Eijndhoven (Utrecht University, The Netherlands).
+//
+// For the minimisation process the TFitter facility, which is basically Minuit,
+// is used. Minimisation is performed by invokation of the SIMPLEX method,
+// followed by an invokation of HESSE to determine the uncertainties on the results. 
+// The statistics of the TFitter result are stored as an AliSignal object
+// in the track, which can be obtained via the GetFitDetails memberfunction.
+// In the minimisation procedure an overall plausibility for the fitted track
+// is determined based on a convoluted Pandel pdf value for each used hit.
+// This track plausibility is expressed in terms of a Bayesian psi value
+// w.r.t. a Convoluted Pandel PDF.
+// The Baysian psi value is defined as -loglikelihood in a decibel scale.
+// This implies psi=-10*log10(L) where L=p(D|HI) being the likelihood of
+// the data D under the hypothesis H and prior information I.
+// Since all (associated) hits contribute independently to the Bayesian psi
+// value, this psi value is built up by summation of the various hit contributions.
+// As such, the FitDetails entries contain the statistics of all the different
+// hit contributions, like PsiMedian, PsiMean, and PsiSigma.
+// The Bayesian psi value is available in the fit details under the name "PsiSum".
+// In addition the standard Minuit results like IERFIT, FCN, EDM etc... are
+// also available from the FitDetails.
+//
+// The convoluted Pandel value is evaluated in various areas in the distance-time
+// space as described in the CPandel writeup by O. Fadiran, G. Japaridze
+// and N. van Eijndhoven.
+// In case the distance-time point of a certain hit falls outside the
+// validity rectangle, the point is moved onto the corresponding side location
+// of the rectangle. For this new location the Pandel value is evaluated for
+// this hit and an extra penalty is added to the corresponding psi value
+// for this hit. 
+// By default this penalty value amounts to 0 dB, but the user can
+// modify this penalty value via the memberfunction SetPenalty.
+// This allows investigation/tuning of the sensitivity to hits with
+// extreme distance and/or time residual values.
 //
 // Use the UseTracks memberfunction to specify the first guess tracks
 // to be processed by the minimiser.
 // Use the SelectHits memberfunction to specify the hits to be used.
 // By default only the hits associated to the first guess track are used.
 //
-// The fit processor (TFitter which is basically Minuit) printlevel
-// can be selected via the memberfunction SetPrintLevel.
+// The fit processor printlevel can be selected via the memberfunction SetPrintLevel.
 // By default all printout is suppressed (i.e. level=-2).
-//
-// In case of bad parameter values as input for the minimiser, the
-// value of the Pandel sometimes cannot be evaluated.
-// In such a case a penalty value will be added.
-// By default this penalty value amounts to 3 dB, but the user can
-// modify this penalty value via the memberfunction SetPenalty.
-// This allows investigation/tuning of the sensitivity to hits with
-// bad time residuals.
 // 
 // An example of how to invoke this processor after Xtalk, hit cleaning
 // and a direct walk first guess estimate can be found in the ROOT example
 // via the GetParentTrack facility of these "IcePandel" tracks.
 // Furthermore, all the hits that were used in the minisation are available
 // via the GetSignal facility of a certain track.
-// The statistics of the TFitter result are stored as an AliSignal object
-// in the track, which can be obtained via the GetFitDetails memberfunction.
-// Currently no overall probability has yet been defined to indicate the
-// quality of a certain fit result. So, for the time being the user has
-// to judge the fit quality his/herself by means of the various TFitter stats
-// as accessible via the GetFitDetails facility.
 //
 // An example of how the various data can be accessed is given below,
 // where "evt" indicates the pointer to the IceEvent structure.
@@ -134,12 +154,13 @@ IcePandel::IcePandel(const char* name,const char* title) : TTask(name,title)
  fEvt=0;
  fUseNames=0;
  fUseNtk=0;
- fTrack=0;
  fHits=0;
  fFitter=0;
  fTrackname="IcePandel";
  fCharge=0;
- fPenalty=3;
+ fPenalty=0;
+ fTkfit=0;
+ fFitstats=0;
 
  // Set the global pointer to this instance
  gIcePandel=this;
@@ -168,6 +189,16 @@ IcePandel::~IcePandel()
   delete fFitter;
   fFitter=0;
  }
+ if (fTkfit)
+ {
+  delete fTkfit;
+  fTkfit=0;
+ }
+ if (fFitstats)
+ {
+  delete fFitstats;
+  fFitstats=0;
+ }
 }
 ///////////////////////////////////////////////////////////////////////////
 void IcePandel::Exec(Option_t* opt)
@@ -203,10 +234,14 @@ void IcePandel::Exec(Option_t* opt)
   cout << " *IcePandel* Hit selection mode : " << fSelhits << endl;
   cout << " *IcePandel* Penalty value for minimiser : " << fPenalty << " dB." << endl;
   cout << endl;
+
+  fPsistats.SetStoreMode(1);
+
   fFirst=0;
  }
 
- Double_t pi=acos(-1.);
+ const Double_t pi=acos(-1.);
+ const Double_t e=exp(1.);
 
  // Initialisation of the minimisation processor
  Double_t arglist[100];
@@ -240,23 +275,33 @@ void IcePandel::Exec(Option_t* opt)
  // Track by track processing of the selected first guess classes
  Float_t vec[3];
  Float_t err[3];
- Float_t margin,minmarg=25; // (minimal) margin for the fitter r0 search area
  Ali3Vector p;
  AliPosition pos;
- AliTrack tkfit;
- tkfit.SetNameTitle(fTrackname.Data(),"IcePandel fit result");
- AliSignal fitstats;
- fitstats.SetNameTitle("Fitstats","TFitter stats for Pandel fit");
- fitstats.SetSlotName("IER",1);
- fitstats.SetSlotName("FCN",2);
- fitstats.SetSlotName("EDM",3);
- fitstats.SetSlotName("NVARS",4);
- Float_t x,y,z,theta,phi;
- Float_t xmin,xmax,ymin,ymax,zmin,zmax,thetamin,thetamax,phimin,phimax;
+ if (!fTkfit)
+ {
+  fTkfit=new AliTrack();
+  fTkfit->SetNameTitle(fTrackname.Data(),"IcePandel fit result");
+ }
+ if (!fFitstats)
+ {
+  fFitstats=new AliSignal();
+  fFitstats->SetNameTitle("Fitstats","TFitter stats for Pandel fit");
+  fFitstats->SetSlotName("IERFIT",1);
+  fFitstats->SetSlotName("FCN",2);
+  fFitstats->SetSlotName("EDM",3);
+  fFitstats->SetSlotName("NVARS",4);
+  fFitstats->SetSlotName("IERERR",5);
+  fFitstats->SetSlotName("PsiSum",6);
+  fFitstats->SetSlotName("PsiMedian",7);
+  fFitstats->SetSlotName("PsiMean",8);
+  fFitstats->SetSlotName("PsiSigma",9);
+ }
+ Float_t x,y,z,theta,phi,t0;
  Double_t amin,edm,errdef; // Minimisation stats
- Int_t ier,nvpar,nparx;    // Minimisation stats
+ Int_t ierfit,iererr,nvpar,nparx;    // Minimisation stats
  Int_t ntk=0;
  Int_t nsig=0;
+ AliTrack* track=0;
  for (Int_t iclass=0; iclass<nclasses; iclass++) // Loop over first guess classes
  {
   strx=(TObjString*)fUseNames->At(iclass);
@@ -269,20 +314,22 @@ void IcePandel::Exec(Option_t* opt)
 
   for (Int_t jtk=0; jtk<ntk; jtk++) // Loop over tracks of a certain class
   {
-   fTrack=(AliTrack*)tracks->At(jtk);
-   if (!fTrack) continue;
+   track=(AliTrack*)tracks->At(jtk);
+   if (!track) continue;
 
-   AliPosition* r0=fTrack->GetReferencePoint();
+   AliPosition* r0=track->GetReferencePoint();
    if (!r0) continue;
 
+   AliTimestamp* tt0=r0->GetTimestamp();
+
    // If selected, use only the first guess track associated hits
    if (fSelhits==1)
    {
     fHits->Clear();
-    nsig=fTrack->GetNsignals();
+    nsig=track->GetNsignals();
     for (Int_t is=1; is<=nsig; is++)
     {
-     AliSignal* sx=fTrack->GetSignal(is);
+     AliSignal* sx=track->GetSignal(is);
      if (!sx) continue;
      if (!sx->GetDevice()->InheritsFrom("IceGOM")) continue;
      if (sx->GetDeadValue("ADC") || sx->GetDeadValue("LE") || sx->GetDeadValue("TOT")) continue;
@@ -298,34 +345,19 @@ void IcePandel::Exec(Option_t* opt)
    x=vec[0];
    y=vec[1];
    z=vec[2];
-   margin=5.*err[0];
-   if (margin<minmarg) margin=minmarg;
-   xmin=x-margin;
-   xmax=x+margin;
-   margin=5.*err[1];
-   if (margin<minmarg) margin=minmarg;
-   ymin=y-margin;
-   ymax=y+margin;
-   margin=5.*err[2];
-   if (margin<minmarg) margin=minmarg;
-   zmin=z-margin;
-   zmax=z+margin;
-
-   p=fTrack->Get3Momentum();
+
+   p=track->Get3Momentum();
    p.GetVector(vec,"sph");
 
    theta=vec[1];
    phi=vec[2];
-   thetamin=theta-(pi/3.);
-   if (thetamin<0) thetamin=0;
-   thetamax=theta+(pi/3.);
-   if (thetamax>pi) thetamax=pi;
-   phimin=phi-(pi/2.);
-   phimax=phi+(pi/2.);
 
+   t0=fEvt->GetDifference(tt0,"ns");
+   
    // Process this first guess track with its associated hits
    fFitter->Clear();
 
+   // Set user selected TFitter printout level
    arglist[0]=fPrint;
    if (fPrint==-2) arglist[0]=-1;
    fFitter->ExecuteCommand("SET PRINT",arglist,1);
@@ -333,25 +365,41 @@ void IcePandel::Exec(Option_t* opt)
 
    fFitter->SetFitMethod("loglikelihood");
 
-   fFitter->SetParameter(0,"r0x",x,0.001,xmin,xmax);
-   fFitter->SetParameter(1,"r0y",y,0.001,ymin,ymax);
-   fFitter->SetParameter(2,"r0z",z,0.001,zmin,zmax);
-   fFitter->SetParameter(3,"theta",theta,0.001,thetamin,thetamax);
-   fFitter->SetParameter(4,"phi",phi,0.001,phimin,phimax);
+   // Define errors to represent 1 sigma for this likelihood scale
+   arglist[0]=5.*log10(e);
+   fFitter->ExecuteCommand("SET ERRORDEF",arglist,1);
+
+   fFitter->SetParameter(0,"r0x",x,0.1,0,0);
+   fFitter->SetParameter(1,"r0y",y,0.1,0,0);
+   fFitter->SetParameter(2,"r0z",z,0.1,0,0);
+   fFitter->SetParameter(3,"theta",theta,0.001,0,pi);
+   fFitter->SetParameter(4,"phi",phi,0.001,0,2.*pi);
+   fFitter->SetParameter(5,"t0",t0,1.,0,32000);
 
    fFitter->SetFCN(IcePandelFCN);
 
+   fTkfit->Reset();
+
    arglist[0]=0;
-   ier=fFitter->ExecuteCommand("MINIMIZE",arglist,0);
+   ierfit=fFitter->ExecuteCommand("SIMPLEX",arglist,0);
 
    fFitter->GetStats(amin,edm,errdef,nvpar,nparx);
-   fitstats.Reset();
-   fitstats.SetSignal(ier,1);
-   fitstats.SetSignal(amin,2);
-   fitstats.SetSignal(edm,3);
-   fitstats.SetSignal(nvpar,4);
 
-   // Resulting parameters after minimisation
+   fFitstats->Reset();
+   fFitstats->SetSignal(ierfit,1);
+   fFitstats->SetSignal(amin,2);
+   fFitstats->SetSignal(edm,3);
+   fFitstats->SetSignal(nvpar,4);
+
+   fFitstats->SetSignal(fPsistats.GetSum(1),6);
+   fFitstats->SetSignal(fPsistats.GetMedian(1),7);
+   fFitstats->SetSignal(fPsistats.GetMean(1),8);
+   fFitstats->SetSignal(fPsistats.GetSigma(1),9);
+
+   iererr=fFitter->ExecuteCommand("HESSE",arglist,0);
+   fFitstats->SetSignal(iererr,5);
+
+   // Resulting parameters after minimisation and error calculation
    vec[0]=fFitter->GetParameter(0);
    vec[1]=fFitter->GetParameter(1);
    vec[2]=fFitter->GetParameter(2);
@@ -370,24 +418,26 @@ void IcePandel::Exec(Option_t* opt)
    p.SetVector(vec,"sph");
    p.SetErrors(err,"sph");
 
+   t0=fFitter->GetParameter(5);
+   AliTimestamp t0fit((AliTimestamp)(*fEvt));
+   t0fit.Add(0,0,t0);
+
    // Enter the fit result as a track in the event structure
-   tkfit.Reset();
    ntkreco++;
-   tkfit.SetId(ntkreco);
-   tkfit.SetCharge(fCharge);
-   tkfit.SetParentTrack(fTrack);
-   AliTimestamp* tt0=r0->GetTimestamp();
-   pos.SetTimestamp(*tt0);
-   tkfit.SetTimestamp(*tt0);
-   tkfit.SetReferencePoint(pos);
-   tkfit.Set3Momentum(p);
-   tkfit.SetFitDetails(fitstats);
+   fTkfit->SetId(ntkreco);
+   fTkfit->SetCharge(fCharge);
+   fTkfit->SetParentTrack(track);
+   pos.SetTimestamp(t0fit);
+   fTkfit->SetTimestamp(t0fit);
+   fTkfit->SetReferencePoint(pos);
+   fTkfit->Set3Momentum(p);
    for (Int_t ihit=0; ihit<fHits->GetEntries(); ihit++)
    {
     AliSignal* sx=(AliSignal*)fHits->At(ihit);
-    if (sx) tkfit.AddSignal(*sx);
+    if (sx) fTkfit->AddSignal(*sx);
    }
-   fEvt->AddTrack(tkfit);
+   fTkfit->SetFitDetails(fFitstats);
+   fEvt->AddTrack(fTkfit);
   } // End loop over tracks
  } // End loop over first guess classes
 
@@ -422,10 +472,10 @@ void IcePandel::UseTracks(TString classname,Int_t n)
 // Example :
 // ---------
 // UseTracks("IceDwalk",5);
-// UseTracks("IceLfit",2);
+// UseTracks("IceLinefit",2);
 // UseTracks("IceJams");
 //
-// This will use the first 5 IceDwalk, the first 2 IceLfit and all the
+// This will use the first 5 IceDwalk, the first 2 IceLinefit and all the
 // IceJams tracks which are encountered in the event structure.
 
  if (!fUseNames)
@@ -490,16 +540,20 @@ void IcePandel::SetCharge(Float_t charge)
 ///////////////////////////////////////////////////////////////////////////
 void IcePandel::SetPenalty(Float_t val)
 {
-// Set user defined penalty value in dB for the minimiser outside the range.
-// This allows investigation/tuning of the sensitivity to hits with bad
-// time residuals.
-// By default the penalty val=3 dB is set in the constructor of this class.
+// Set user defined psi penalty value (in dB) in the minimiser for
+// distance-time points that fall outside the validity rectangle.
+// This allows investigation/tuning of the sensitivity to hits with
+// extreme distance and/or time residual values.
+// By default the penalty val=0 is set in the constructor of this class.
  fPenalty=val;
 }
 ///////////////////////////////////////////////////////////////////////////
 void IcePandel::FitFCN(Int_t&,Double_t*,Double_t& f,Double_t* x,Int_t)
 {
-// The Pandel function used for the minimisation process.
+// Minimisation of the Bayesian psi value for a track w.r.t. a Convoluted Pandel PDF.
+// The Baysian psi value is defined as -loglikelihood in a decibel scale.
+// This implies psi=-10*log10(L) where L=p(D|HI) being the likelihood of
+// the data D under the hypothesis H and prior information I.
 
  const Float_t c=0.299792;           // Light speed in vacuum in meters per ns
  const Float_t nice=1.35634;         // Refractive index of ice
@@ -509,17 +563,14 @@ void IcePandel::FitFCN(Int_t&,Double_t*,Double_t& f,Double_t* x,Int_t)
  const Float_t cice=c/nice;          // Light speed in ice in meters per ns
  const Float_t tau=557;
  const Double_t rho=((1./tau)+(cice/labs));
- const Double_t e=exp(1.);
+ const Double_t pi=acos(-1.);
 
- f=0;
+ // Assumed PMT timing jitter in ns
+ const Double_t sigma=10;
 
- // The first original guess track parameters
- AliPosition* tr0=fTrack->GetReferencePoint();
- AliTimestamp* tt0=tr0->GetTimestamp();
- Float_t t0=fEvt->GetDifference(tt0,"ns");
- Ali3Vector tp=fTrack->Get3Momentum();
+ f=0;
 
- // The new r0 and p vectors from the minimisation
+ // The new r0 and p vectors and t0 from the minimisation
  Float_t vec[3];
  
  AliPosition r0;
@@ -534,40 +585,132 @@ void IcePandel::FitFCN(Int_t&,Double_t*,Double_t& f,Double_t* x,Int_t)
  vec[2]=x[4];
  p.SetVector(vec,"sph");
 
+ Float_t t0=x[5];
+
+ // Construct a track with the new values from the minimisation
+ fTkfit->SetReferencePoint(r0);
+ fTkfit->Set3Momentum(p);
+
  Int_t nhits=fHits->GetEntries();
+
  AliPosition rhit;
  Ali3Vector r12;
  Float_t d,dist,thit,tgeo;
- Double_t tres,zeta,pandel;
- for (Int_t i=0; i<nhits; i++)
+ Double_t tres,ksi,eta,pandel;
+ Double_t cpandel1,cpandel2,cpandel3,cpandel;
+ Double_t z,k,alpha,beta,phi,n1,n2,n3,u; // Function parameters for regions 3 and 4
+ Int_t ier;
+ Double_t psihit=0;
+ fPsistats.Reset();
+ for (Int_t i=1; i<=nhits; i++)
  {
   AliSignal* sx=(AliSignal*)fHits->At(i);
   if (!sx) continue;
   IceGOM* omx=(IceGOM*)sx->GetDevice();
   if (!omx) continue;
   rhit=omx->GetPosition();
-  d=fTrack->GetDistance(rhit);
-  zeta=d/lambda;
-  d*=sin(thetac);
+  d=fTkfit->GetDistance(rhit);
+  ksi=d/lambda;
   r12=rhit-r0;
   dist=p.Dot(r12)+d*tan(thetac);
   tgeo=t0+dist/c;
   thit=sx->GetSignal("LE",7);
   tres=thit-tgeo;
 
-  // LLH optimisation based on a decibel scaled Pandel function
-  // Avoid minimiser problem for infinite derivative on Pandel surface
-  // This problem will be absent when using a smooth convoluted Pandel
-  // Use -10*log10(p) expression to obtain intuitive decibel scale for return value "f"
-  if (tres>=0.001 && tres<=10000 && zeta>=0.001 && zeta<=10000)
+  // The Convoluted Pandel function evaluation
+  // For definitions of functions and validity regions, see the
+  // CPandel writeup of O. Fadiran, G. Japaridze and N. van Eijndhoven
+
+  // Move points which are outside the validity rectangle in the (tres,ksi) space
+  // to the edge of the validity rectangle and signal the use of the penalty
+  ier=0;
+  if (tres<-25.*sigma)
+  {
+   tres=-25.*sigma;
+   ier=1;
+  }
+  if (tres>3500)
+  {
+   tres=3500;
+   ier=1;
+  }
+  if (ksi>50)
+  {
+   ksi=50;
+   ier=1;
+  }
+
+  eta=(rho*sigma)-(tres/sigma);
+
+  if (ksi<=0) // The zero distance (ksi=0) axis
+  {
+   cpandel=exp(-tres*tres/(2.*sigma*sigma))/(sigma*sqrt(2.*pi));
+  }
+  else if (ksi<=5 && tres>=-5.*sigma && tres<=30.*sigma) // The exact expression in region 1
+  {
+   cpandel1=pow(rho,ksi)*pow(sigma,ksi-1.)*exp(-tres*tres/(2.*sigma*sigma))/pow(2.,0.5*(1.+ksi));
+   cpandel2=ROOT::Math::conf_hyperg(ksi/2.,0.5,eta*eta/2.)/TMath::Gamma((ksi+1.)/2.);
+   cpandel3=sqrt(2.)*eta*ROOT::Math::conf_hyperg((ksi+1.)/2.,1.5,eta*eta/2.)/TMath::Gamma(ksi/2.);
+
+   cpandel=cpandel1*(cpandel2-cpandel3);
+  }
+  else if (ksi<=1 && tres>30.*sigma && tres<=3500) // Approximation in region 2
+  {
+   pandel=pow(rho,ksi)*pow(tres,(ksi-1.))*exp(-rho*tres)/TMath::Gamma(ksi);
+
+   cpandel=exp(rho*rho*sigma*sigma/2.)*pandel;
+  }
+  else if (ksi<=1 && tres<-5.*sigma && tres>=-25.*sigma) // Approximation in region 5
+  {
+   cpandel=pow(rho*sigma,ksi)*pow(eta,-ksi)*exp(-tres*tres/(2.*sigma*sigma))/(sigma*sqrt(2.*pi));
+  }
+  else if (ksi<=50 && tres>=0 && tres<=3500) // Approximation in region 3
+  {
+   z=-eta/sqrt(4.*ksi-2.);
+   k=0.5*(z*sqrt(1.+z*z)+log(z+sqrt(1.+z*z)));
+   alpha=-tres*tres/(2.*sigma*sigma)+eta*eta/4.-ksi/2.+0.25+k*(2.*ksi-1.);
+   alpha+=-log(1.+z*z)/4.-ksi*log(2.)/2.+(ksi-1.)*log(2.*ksi-1.)/2.+ksi*log(rho)+(ksi-1.)*log(sigma);
+   beta=0.5*(z/sqrt(1.+z*z)-1.);
+   n1=beta*(20.*beta*beta+30.*beta+9.)/12.;
+   n2=pow(beta,2.)*(6160.*pow(beta,4.)+18480.*pow(beta,3.)+19404.*pow(beta,2.)+8028.*beta+945.)/288.;
+   n3=27227200.*pow(beta,6.)+122522400.*pow(beta,5.)+220540320.*pow(beta,4.);
+   n3+=200166120.*pow(beta,3.)+94064328.*pow(beta,2.)+20546550.*beta+1403325.;
+   n3*=pow(beta,3.)/51840.;
+   phi=1.-n1/(2.*ksi-1.)+n2/pow(2.*ksi-1.,2.)-n3/pow(2.*ksi-1.,3.);
+   cpandel=exp(alpha)*phi/TMath::Gamma(ksi);
+  }
+  else if (ksi<=50 && tres<0 && tres>=-25.*sigma) // Approximation in region 4
   {
-   pandel=-10.*((zeta*log10(rho))+((zeta-1.)*log10(tres))-(rho*tres*log10(e))-log10(TMath::Gamma(zeta)));
-   f+=pandel;
+   z=eta/sqrt(4.*ksi-2.);
+   k=0.5*(z*sqrt(1.+z*z)+log(z+sqrt(1.+z*z)));
+   u=exp(ksi/2.-0.25)*pow(2.*ksi-1.,-ksi/2.)*pow(2.,(ksi-1.)/2.);
+   beta=0.5*(z/sqrt(1.+z*z)-1.);
+   n1=beta*(20.*beta*beta+30.*beta+9.)/12.;
+   n2=pow(beta,2.)*(6160.*pow(beta,4.)+18480.*pow(beta,3.)+19404.*pow(beta,2.)+8028.*beta+945.)/288.;
+   n3=27227200.*pow(beta,6.)+122522400.*pow(beta,5.)+220540320.*pow(beta,4.);
+   n3+=200166120.*pow(beta,3.)+94064328.*pow(beta,2.)+20546550.*beta+1403325.;
+   n3*=pow(beta,3.)/51840.;
+   phi=1.+n1/(2.*ksi-1.)+n2/pow(2.*ksi-1.,2.)+n3/pow(2.*ksi-1.,3.);
+   cpandel=pow(rho,ksi)*pow(sigma,ksi-1.)*exp(-pow(tres,2.)/(2.*pow(sigma,2.))+pow(eta,2.)/4.)/sqrt(2.*pi);
+   cpandel*=u*phi*exp(-k*(2.*ksi-1.))*pow(1.+z*z,-0.25);
   }
-  else
+  else // (tres,ksi) outside validity rectangle
   {
-   f+=fPenalty; // Penalty in case tres and/or zeta outside range
+   ier=1;
+   cout << " *IcePandel::FitFCN* Outside rectangle. We should never get here." << endl;
   }
+
+  // Use 10*log10 expression to obtain intuitive dB scale
+  // Omit (small) negative values which are possible due to computer accuracy
+  psihit=0;
+  if (cpandel>0) psihit=-10.*log10(cpandel);
+
+  // Penalty in dB for (tres,ksi) points outside the validity rectangle
+  if (ier) psihit+=fPenalty; 
+
+  // Update the psi statistics for this hit
+  fPsistats.Enter(float(psihit));
+  f+=psihit;
  }
 }
 ///////////////////////////////////////////////////////////////////////////
index 314a79f..35360d1 100644 (file)
 #include "TString.h"
 #include "TObjString.h"
 #include "TFitter.h"
+#include "Math/SpecFuncMathMore.h"
 
 #include "AliJob.h"
+#include "AliSample.h"
 #include "IceEvent.h"
 #include "IceGOM.h"
 
@@ -27,7 +29,7 @@ class IcePandel : public TTask
   void SelectHits(Int_t mode=1);                       // Specify which hits to be used
   void SetTrackName(TString s);  // Set (alternative) name for the produced tracks
   void SetCharge(Float_t charge);// Set user defined charge for the produced tracks
-  void SetPenalty(Float_t val);  // Set penalty value in dB for the minimiser outside the allowed range
+  void SetPenalty(Float_t val);  // Set psi penalty value in dB for extreme distance and/or time values
   void FitFCN(Int_t&,Double_t*,Double_t&,Double_t*,Int_t); // The minimisation FCN
 
  protected :
@@ -37,13 +39,15 @@ class IcePandel : public TTask
   IceEvent* fEvt;       // Pointer to the current event structure
   TObjArray* fUseNames; // The first guess classnames to be used 
   TArrayI* fUseNtk;     // The max. numbers of the various first guess tracks to be used
-  AliTrack* fTrack;     // Pointer to the first guess track being processed
   TObjArray* fHits;     // The various hits to be used in the fitting process 
   TFitter* fFitter;     // Pointer to the minimisation processor
   TString fTrackname;   // The name identifier for the produced tracks
   Float_t fCharge;      // User defined charge of the produced tracks
-  Float_t fPenalty;     // User defined penalty value in dB for the minimiser outside the allowed range
+  Float_t fPenalty;     // User defined psi penalty value (dB) for extreme distance and/or time values
+  AliTrack* fTkfit;     // Pointer to the produced fitted track
+  AliSignal* fFitstats; // The fit details of the produced fitted track
+  AliSample fPsistats;  // Statistics of the Bayesian psi value for the best fitted track
 
- ClassDef(IcePandel,4) // TTask derived class to perform Pandel fitting
+ ClassDef(IcePandel,5) // TTask derived class to perform Convoluted Pandel fitting
 };
 #endif
index 859b86c..2cfe7b1 100644 (file)
@@ -68,3 +68,6 @@
 10-may-2006 NvE Compiler warnings activated for overloaded virtuals in installation
                 script for AMD Athlon in the /scripts sub-directory.
 23-jun-2006 NvE New lightweight AliSignal::AddTrack facility used in IceLinefit and IceDwalk.
+07-jul-2006 NvE Track Id correctly set in IceLinefit.
+28-jul-2006 NvE New class IceChi2 introduced and example macro icechi2.cc added.
+                Also IcePandel updated to use a Convoluted Pandel pdf.
diff --git a/RALICE/icepack/macros/icechi2.cc b/RALICE/icepack/macros/icechi2.cc
new file mode 100644 (file)
index 0000000..f28c521
--- /dev/null
@@ -0,0 +1,70 @@
+//////////////////////////////////////////////////////////
+// Example macro to demonstrate the usage of the IceChi2
+// processor for chi-squared fit reconstruction as a subtask
+// of an F2K conversion job. The macro also shows how one
+// can interactively invoke one or more subtasks
+// (i.e. IceXtalk and EvtAna.cxx) to be executed.
+// The latter is very convenient in testing/comparing
+// new reconstruction/analysis algorithms.
+//
+// To run this macro in batch, just do
+//
+// root -b -q icechi2.cc
+//
+// For more details see the docs of class IceChi2
+//
+// NvE 28-jul-2006 Utrecht University
+//////////////////////////////////////////////////////////
+{
+ gSystem->Load("ralice");
+ gSystem->Load("icepack");
+ gSystem->Load("iceconvert");
+
+ // Interactively compile and load the EvtAna.cxx code
+ gROOT->LoadMacro("EvtAna.cxx+");
+
+ // The database conversion job
+ IceCal2Root cal("IceCal2Root","Calibration format conversion");
+ cal.SetAmacalibFile("amacalib_amanda2_2003.txt");
+ cal.SetOutputFile("cal2003.root");
+ cal.ExecuteJob();
+
+ // The chi-squared fitting procedure
+ IceChi2 chi2("IceChi2","Chi-squared fitting");
+ chi2.UseTracks("IceDwalk");
+
+ // The direct walk reconstruction task
+ IceDwalk dwalk("IceDwalk","Direct walk reconstruction");
+
+ // The calibration processor task
+ IceCalibrate calib("IceCalibrate","Signal calibration");
+ calib.SetCalibFile("cal2003.root");
+
+ // The Xtalk correction processor task
+ IceXtalk xtalk("IceXtalk","Cross talk correction");
+ xtalk.SetCalibFile("cal2003.root");
+
+ // The hit cleaning processor task
+ IceCleanHits clean("IceCleanHits","Hit cleaning");
+
+ // The event analysis task
+ EvtAna evtana("evtana","Event analysis");
+
+ // The F2K event data processing job
+ IceF2k q("IceF2k","Processing of the F2K event data");
+ q.SetMaxEvents(2);
+ q.SetPrintFreq(0);
+ q.SetInputFile("real-reco.f2k");
+ q.SetOutputFile("real-reco.root");
+
+ // Add the various processors as subtasks to the F2K job
+ q.Add(&calib);
+ q.Add(&xtalk);
+ q.Add(&clean);
+ q.Add(&dwalk);
+ q.Add(&chi2);
+ q.Add(&evtana);
+
+ // Perform the conversion
+ q.ExecuteJob();
+}