From: nick Date: Fri, 28 Jul 2006 13:35:54 +0000 (+0000) Subject: 21-jul-2006 NvE AliSample extended for computation of median. X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=219e9b7e7b8c47b2eb80861d12a62bbad0dd613a 21-jul-2006 NvE AliSample extended for computation of median. 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/AliSample.cxx b/RALICE/AliSample.cxx index 62aa9c41e06..5a1abc89068 100644 --- a/RALICE/AliSample.cxx +++ b/RALICE/AliSample.cxx @@ -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; iSet(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; iGetSize() < fN) fArr->Set(fN); + } + + // Order the values of the specified variable + Float_t val=0; + Int_t iadd=0; + for (Int_t j=0; jAt(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=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; +} +/////////////////////////////////////////////////////////////////////////// diff --git a/RALICE/AliSample.h b/RALICE/AliSample.h index d9c960fee2d..9f1ca69984b 100644 --- a/RALICE/AliSample.h +++ b/RALICE/AliSample.h @@ -8,6 +8,7 @@ #include #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. }; diff --git a/RALICE/history.txt b/RALICE/history.txt index 23a58812784..7e229bd8be6 100644 --- a/RALICE/history.txt +++ b/RALICE/history.txt @@ -697,3 +697,4 @@ 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. diff --git a/RALICE/icepack/ICEHeaders.h b/RALICE/icepack/ICEHeaders.h index 2d95d1ef1ac..ac5a431bcaa 100644 --- a/RALICE/icepack/ICEHeaders.h +++ b/RALICE/icepack/ICEHeaders.h @@ -20,4 +20,5 @@ #include "IceDwalk.h" #include "IcePandel.h" #include "IceLinefit.h" +#include "IceChi2.h" diff --git a/RALICE/icepack/ICELinkDef.h b/RALICE/icepack/ICELinkDef.h index 35d99d4a791..79f333ab2e1 100644 --- a/RALICE/icepack/ICELinkDef.h +++ b/RALICE/icepack/ICELinkDef.h @@ -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 index 00000000000..b62ee6339b2 --- /dev/null +++ b/RALICE/icepack/IceChi2.cxx @@ -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; jtkGetEntries(); 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; iAt(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; ihGetEntries(); 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; iclassAt(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; jtkAt(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; ihitGetEntries(); 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; iAt(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; iAt(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 index 00000000000..c02766acbb5 --- /dev/null +++ b/RALICE/icepack/IceChi2.h @@ -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 diff --git a/RALICE/icepack/IceLinefit.cxx b/RALICE/icepack/IceLinefit.cxx index 6ee1191228f..dd38447d88f 100644 --- a/RALICE/icepack/IceLinefit.cxx +++ b/RALICE/icepack/IceLinefit.cxx @@ -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"); diff --git a/RALICE/icepack/IcePandel.cxx b/RALICE/icepack/IcePandel.cxx index 19f53cee580..9f0699c258e 100644 --- a/RALICE/icepack/IcePandel.cxx +++ b/RALICE/icepack/IcePandel.cxx @@ -17,10 +17,45 @@ /////////////////////////////////////////////////////////////////////////// // 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. @@ -29,17 +64,8 @@ // 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 @@ -63,12 +89,6 @@ // 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; iclassAt(iclass); @@ -269,20 +314,22 @@ void IcePandel::Exec(Option_t* opt) for (Int_t jtk=0; jtkAt(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 (marginGet3Momentum(); + + 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; ihitGetEntries(); 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; iAt(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; } } /////////////////////////////////////////////////////////////////////////// diff --git a/RALICE/icepack/IcePandel.h b/RALICE/icepack/IcePandel.h index 314a79f3198..35360d176cc 100644 --- a/RALICE/icepack/IcePandel.h +++ b/RALICE/icepack/IcePandel.h @@ -11,8 +11,10 @@ #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 diff --git a/RALICE/icepack/history.txt b/RALICE/icepack/history.txt index 859b86cca55..2cfe7b1da8f 100644 --- a/RALICE/icepack/history.txt +++ b/RALICE/icepack/history.txt @@ -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 index 00000000000..f28c5216eea --- /dev/null +++ b/RALICE/icepack/macros/icechi2.cc @@ -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(); +}