]> git.uio.no Git - u/mrichter/AliRoot.git/blob - RALICE/icepack/IceChi2.cxx
25-sep-2006 NvE AliSample extended with memberfunction GetSpread() and also
[u/mrichter/AliRoot.git] / RALICE / icepack / IceChi2.cxx
1 /*******************************************************************************
2  * Copyright(c) 2003, IceCube Experiment at the South Pole. All rights reserved.
3  *
4  * Author: The IceCube RALICE-based Offline Project.
5  * Contributors are mentioned in the code where appropriate.
6  *
7  * Permission to use, copy, modify and distribute this software and its
8  * documentation strictly for non-commercial purposes is hereby granted
9  * without fee, provided that the above copyright notice appears in all
10  * copies and that both the copyright notice and this permission notice
11  * appear in the supporting documentation.
12  * The authors make no claims about the suitability of this software for
13  * any purpose. It is provided "as is" without express or implied warranty.
14  *******************************************************************************/
15
16 // $Id$
17
18 ///////////////////////////////////////////////////////////////////////////
19 // Class IceChi2
20 // TTask derived class to perform track fitting via chi-squared minimisation.
21 //
22 // For the minimisation process the TFitter facility, which is basically Minuit,
23 // is used. Minimisation is performed by invokation of the SIMPLEX method,
24 // followed by an invokation of HESSE to determine the uncertainties on the results.
25 // The statistics of the TFitter result are stored as an AliSignal object
26 // in the track, which can be obtained via the GetFitDetails memberfunction.
27 // After the chi-squared minimisation procedure has been performed, an overall
28 // plausibility for the fitted track will be determined based on a convoluted
29 // Pandel pdf value for each used hit.
30 // This track plausibility is expressed in terms of a Bayesian psi value
31 // w.r.t. a Convoluted Pandel PDF.
32 // The Baysian psi value is defined as -loglikelihood in a decibel scale.
33 // This implies psi=-10*log10(L) where L=p(D|HI) being the likelihood of
34 // the data D under the hypothesis H and prior information I.
35 // Since all (associated) hits contribute independently to the Bayesian psi
36 // value, this psi value is built up by summation of the various hit contributions.
37 // As such, the FitDetails entries contain the statistics of all the different
38 // hit contributions, like PsiMedian, PsiMean, and PsiSigma.
39 // The Bayesian psi value is available in the fit details under the name "PsiSum".
40 // In addition the standard Minuit results like IERFIT, FCN, EDM etc... are
41 // also available from the FitDetails.
42 //
43 // The convoluted Pandel value is evaluated in various areas in the distance-time
44 // space as described in the CPandel writeup by O. Fadiran, G. Japaridze
45 // and N. van Eijndhoven.
46 // In case the distance-time point of a certain hit falls outside the
47 // validity rectangle, the point is moved onto the corresponding side location
48 // of the rectangle. For this new location the Pandel value is evaluated for
49 // this hit and an extra penalty is added to the corresponding psi value
50 // for this hit. 
51 // By default this penalty value amounts to 0 dB, but the user can
52 // modify this penalty value via the memberfunction SetPenalty.
53 // This allows investigation/tuning of the sensitivity to hits with
54 // extreme distance and/or time residual values.
55 //
56 // A separate treatment of the phase and group velocities is introduced
57 // which will provide more accurate time residuals due to the different
58 // velocities of the Cerenkov wave front (v_phase) and the actually detected
59 // photons (v_group).
60 // This distinction between v_phase and v_group can be (de)activated via the
61 // memberfunction SetVgroupUsage(). By default the distinction between v_phase
62 // and v_group is activated in the constructor of this class.
63 //
64 // Use the UseTracks memberfunction to specify the first guess tracks
65 // to be processed by the minimiser.
66 // By default only the first encountered IceDwalk track will be processed.
67 //
68 // Use the SelectHits memberfunction to specify the hits to be used.
69 // By default only the hits associated to the first guess track are used.
70 //
71 // Information about the actual parameter settings can be found in the event
72 // structure itself via the device named "IceChi2".
73 //
74 // The fit processor printlevel can be selected via the memberfunction
75 // SetPrintLevel. By default all printout is suppressed (i.e. level=-2).
76 //
77 // An example of how to invoke this processor after Xtalk, hit cleaning
78 // and a direct walk first guess estimate can be found in the ROOT example
79 // macro icechi2.cc which resides in the /macros subdirectory.
80 // 
81 // The minimisation results are stored in the IceEvent structure as
82 // tracks with as default the name "IceChi2" (just like the first guess
83 // results of e.g. IceDwalk).
84 // This track name identifier can be modified by the user via the
85 // SetTrackName() memberfunction. This will allow unique identification
86 // of tracks which are produced when re-processing existing data with
87 // different criteria.
88 // By default the charge of the produced tracks is set to 0, since
89 // no distinction can be made between positive or negative tracks.
90 // However, the user can define the track charge by invokation
91 // of the memberfunction SetCharge().
92 // This facility may be used to distinguish tracks produced by the
93 // various reconstruction algorithms in a (3D) colour display
94 // (see the class AliHelix for further details).  
95 // A pointer to the first guess track which was used as input is available
96 // via the GetParentTrack facility of these "IceChi2" tracks.
97 // Furthermore, all the hits that were used in the minisation are available
98 // via the GetSignal facility of a certain track.
99 //
100 // An example of how the various data can be accessed is given below,
101 // where "evt" indicates the pointer to the IceEvent structure.
102 //
103 // Example for accessing data :
104 // ----------------------------
105 // TObjArray* tracks=evt->GetTracks("IceChi2");
106 // if (!tracks) return;
107 // AliPosition* r0=0;
108 // Float_t psi=0;
109 // for (Int_t jtk=0; jtk<tracks->GetEntries(); jtk++)
110 // {
111 //  AliTrack* tx=(AliTrack*)tracks->At(jtk);
112 //  if (!tx) continue;
113 //  tx->Data("sph");
114 //  r0=tx->GetReferencePoint();
115 //  if (r0) r0->Data();
116 //  sx=(AliSignal*)tx->GetFitDetails();
117 //  if (sx) psi=sx->GetSignal("PsiSum");
118 //  AliTrack* tx2=tx->GetParentTrack();
119 //  if (!tx2) continue;
120 //  tx2->Data("sph");
121 //  r0=tx2->GetReferencePoint();
122 //  if (r0) r0->Data();
123 // }
124 //
125 // Notes :
126 // -------
127 // 1) This processor only works properly on data which are Time and ADC
128 //    calibrated and contain tracks from first guess algorithms like
129 //    e.g. IceDwalk.
130 // 2) In view of the usage of TFitter/Minuit minimisation, a global pointer
131 //    to the instance of this class (gIceChi2) and a global static
132 //    wrapper function (IceChi2FCN) have been introduced, to allow the
133 //    actual minimisation to be performed via the memberfunction FitFCN.
134 //    This implies that in a certain processing job only 1 instance of
135 //    this IceChi2 class may occur.
136 //
137 //--- Author: Nick van Eijndhoven 16-may-2006 Utrecht University
138 //- Modified: NvE $Date$ Utrecht University
139 ///////////////////////////////////////////////////////////////////////////
140  
141 #include "IceChi2.h"
142 #include "Riostream.h"
143
144 // Global pointer to the instance of this object
145  IceChi2* gIceChi2=0;
146
147 // TFitter/Minuit interface to IceChi2::FitFCN
148  void IceChi2FCN(Int_t& npar,Double_t* gin,Double_t& f,Double_t* u,Int_t flag)
149  {
150   if (gIceChi2) gIceChi2->FitFCN(npar,gin,f,u,flag);
151  }
152
153 ClassImp(IceChi2) // Class implementation to enable ROOT I/O
154
155 IceChi2::IceChi2(const char* name,const char* title) : TTask(name,title)
156 {
157 // Default constructor.
158  fFirst=1;
159  fPrint=-2;
160  fSelhits=1;
161  fVgroup=1;
162  fEvt=0;
163  fUseNames=0;
164  fUseNtk=0;
165  fHits=0;
166  fFitter=0;
167  fTrackname="IceChi2";
168  fCharge=0;
169  fPenalty=0;
170  fTkfit=0;
171  fFitstats=0;
172
173  // Set the global pointer to this instance
174  gIceChi2=this;
175 }
176 ///////////////////////////////////////////////////////////////////////////
177 IceChi2::~IceChi2()
178 {
179 // Default destructor.
180  if (fUseNames)
181  {
182   delete fUseNames;
183   fUseNames=0;
184  }
185  if (fUseNtk)
186  {
187   delete fUseNtk;
188   fUseNtk=0;
189  }
190  if (fHits)
191  {
192   delete fHits;
193   fHits=0;
194  }
195  if (fFitter)
196  {
197   delete fFitter;
198   fFitter=0;
199  }
200  if (fTkfit)
201  {
202   delete fTkfit;
203   fTkfit=0;
204  }
205  if (fFitstats)
206  {
207   delete fFitstats;
208   fFitstats=0;
209  }
210 }
211 ///////////////////////////////////////////////////////////////////////////
212 void IceChi2::Exec(Option_t* opt)
213 {
214 // Implementation of the hit fitting procedure.
215
216  TString name=opt;
217  AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
218
219  if (!parent) return;
220
221  fEvt=(IceEvent*)parent->GetObject("IceEvent");
222  if (!fEvt) return;
223
224  // Storage of the used parameters in the IceChi2 device
225  AliSignal params;
226  params.SetNameTitle("IceChi2","IceChi2 processor parameters");
227  params.SetSlotName("Selhits",1);
228  params.SetSlotName("Penalty",2);
229  params.SetSlotName("Vgroup",3);
230
231  params.SetSignal(fSelhits,1);
232  params.SetSignal(fPenalty,2);
233  params.SetSignal(fVgroup,3);
234
235  fEvt->AddDevice(params);
236
237  if (!fUseNames) UseTracks("IceDwalk",1);
238
239  Int_t nclasses=fUseNames->GetEntries(); // Number of first guess classes to be processed
240  Int_t ntkmax=0; // Max. number of tracks for a certain class
241  TObjString* strx=0;
242  TString str;
243
244  if (fFirst)
245  {
246   cout << " *IceChi2* First guess selections to be processed (-1=all)." << endl;
247   for (Int_t i=0; i<nclasses; i++)
248   {
249    strx=(TObjString*)fUseNames->At(i);
250    if (!strx) continue;
251    str=strx->GetString();
252    ntkmax=fUseNtk->At(i);
253    cout << " Maximally " << ntkmax << " track(s) per event for procedure : " << str.Data() << endl;
254   }
255   cout << " *IceChi2* Hit selection mode : " << fSelhits << endl;
256   cout << " *IceChi2* Penalty value for psi evaluation outside range : " << fPenalty << endl;
257   cout << endl;
258
259   fPsistats.SetStoreMode(1);
260
261   fFirst=0;
262  }
263
264  const Double_t pi=acos(-1.);
265
266  // Initialisation of the minimisation processor
267  Double_t arglist[100];
268  if (!fFitter) fFitter=new TFitter();
269
270  // The number of reconstructed tracks already present in the event
271  Int_t ntkreco=fEvt->GetNtracks(1);
272
273  if (!fHits)
274  {
275   fHits=new TObjArray();
276  }
277  else
278  {
279   fHits->Clear();
280  }
281
282  // If selected, use all the good quality hits of the complete event
283  if (fSelhits==0)
284  {
285   TObjArray* hits=fEvt->GetHits("IceGOM");
286   for (Int_t ih=0; ih<hits->GetEntries(); ih++)
287   {
288    AliSignal* sx=(AliSignal*)hits->At(ih);
289    if (!sx) continue;
290    if (sx->GetDeadValue("ADC") || sx->GetDeadValue("LE") || sx->GetDeadValue("TOT")) continue;
291    fHits->Add(sx);
292   }
293  }
294
295  // Track by track processing of the selected first guess classes
296  Float_t vec[3];
297  Float_t err[3];
298  Ali3Vector p;
299  AliPosition pos;
300  if (!fTkfit)
301  {
302   fTkfit=new AliTrack();
303   fTkfit->SetNameTitle(fTrackname.Data(),"IceChi2 fit result");
304  }
305  if (!fFitstats)
306  {
307   fFitstats=new AliSignal();
308   fFitstats->SetNameTitle("Fitstats","TFitter stats for Chi2 fit");
309   fFitstats->SetSlotName("IERFIT",1);
310   fFitstats->SetSlotName("FCN",2);
311   fFitstats->SetSlotName("EDM",3);
312   fFitstats->SetSlotName("NVARS",4);
313   fFitstats->SetSlotName("IERERR",5);
314   fFitstats->SetSlotName("PsiSum",6);
315   fFitstats->SetSlotName("PsiMedian",7);
316   fFitstats->SetSlotName("PsiSpread",8);
317   fFitstats->SetSlotName("PsiMean",9);
318   fFitstats->SetSlotName("PsiSigma",10);
319  }
320  Float_t x,y,z,theta,phi,t0;
321  Double_t amin,edm,errdef; // Minimisation stats
322  Int_t ierfit,iererr,nvpar,nparx; // Minimisation stats
323  Double_t psi; // Bayesian psi value for the fitted track w.r.t. a Convoluted Pandel PDF
324  Int_t ntk=0;
325  Int_t nsig=0;
326  AliTrack* track=0;
327  for (Int_t iclass=0; iclass<nclasses; iclass++) // Loop over first guess classes
328  {
329   strx=(TObjString*)fUseNames->At(iclass);
330   if (!strx) continue;
331   str=strx->GetString();
332   ntkmax=fUseNtk->At(iclass);
333   TObjArray* tracks=fEvt->GetTracks(str);
334   ntk=tracks->GetEntries();
335   if (ntkmax>0 && ntk>ntkmax) ntk=ntkmax;
336
337   for (Int_t jtk=0; jtk<ntk; jtk++) // Loop over tracks of a certain class
338   {
339    track=(AliTrack*)tracks->At(jtk);
340    if (!track) continue;
341
342    AliPosition* r0=track->GetReferencePoint();
343    if (!r0) continue;
344
345    AliTimestamp* tt0=r0->GetTimestamp();
346
347    // If selected, use only the first guess track associated hits
348    if (fSelhits==1)
349    {
350     fHits->Clear();
351     nsig=track->GetNsignals();
352     for (Int_t is=1; is<=nsig; is++)
353     {
354      AliSignal* sx=track->GetSignal(is);
355      if (!sx) continue;
356      if (!sx->GetDevice()->InheritsFrom("IceGOM")) continue;
357      if (sx->GetDeadValue("ADC") || sx->GetDeadValue("LE") || sx->GetDeadValue("TOT")) continue;
358      fHits->Add(sx);
359     }
360    }
361    
362    if (!fHits->GetEntries()) continue;
363
364    r0->GetVector(vec,"car");
365    r0->GetErrors(err,"car");
366
367    x=vec[0];
368    y=vec[1];
369    z=vec[2];
370
371    p=track->Get3Momentum();
372    p.GetVector(vec,"sph");
373
374    theta=vec[1];
375    phi=vec[2];
376
377    t0=fEvt->GetDifference(tt0,"ns");
378    
379    // Process this first guess track with its associated hits
380    fFitter->Clear();
381
382    // Set user selected TFitter printout level
383    arglist[0]=fPrint;
384    if (fPrint==-2) arglist[0]=-1;
385    fFitter->ExecuteCommand("SET PRINT",arglist,1);
386    if (fPrint==-2) fFitter->ExecuteCommand("SET NOWARNINGS",arglist,0);
387
388    fFitter->SetFitMethod("chisquare");
389
390    fFitter->SetParameter(0,"r0x",x,0.1,0,0);
391    fFitter->SetParameter(1,"r0y",y,0.1,0,0);
392    fFitter->SetParameter(2,"r0z",z,0.1,0,0);
393    fFitter->SetParameter(3,"theta",theta,0.001,0,pi);
394    fFitter->SetParameter(4,"phi",phi,0.001,0,2.*pi);
395    fFitter->SetParameter(5,"t0",t0,1.,0,32000);
396
397    fFitter->SetFCN(IceChi2FCN);
398
399    fTkfit->Reset();
400
401    arglist[0]=0;
402    ierfit=fFitter->ExecuteCommand("SIMPLEX",arglist,0);
403
404    fFitter->GetStats(amin,edm,errdef,nvpar,nparx);
405    fFitstats->Reset();
406    fFitstats->SetSignal(ierfit,1);
407    fFitstats->SetSignal(amin,2);
408    fFitstats->SetSignal(edm,3);
409    fFitstats->SetSignal(nvpar,4);
410
411    iererr=fFitter->ExecuteCommand("HESSE",arglist,0);
412    fFitstats->SetSignal(iererr,5);
413
414
415    // Resulting parameters after minimisation and error calculation
416    vec[0]=fFitter->GetParameter(0);
417    vec[1]=fFitter->GetParameter(1);
418    vec[2]=fFitter->GetParameter(2);
419    err[0]=fFitter->GetParError(0);
420    err[1]=fFitter->GetParError(1);
421    err[2]=fFitter->GetParError(2);
422    pos.SetPosition(vec,"car");
423    pos.SetPositionErrors(err,"car");
424
425    vec[0]=1.;
426    vec[1]=fFitter->GetParameter(3);
427    vec[2]=fFitter->GetParameter(4);
428    err[0]=0.;
429    err[1]=fFitter->GetParError(3);
430    err[2]=fFitter->GetParError(4);
431    p.SetVector(vec,"sph");
432    p.SetErrors(err,"sph");
433
434    t0=fFitter->GetParameter(5);
435    AliTimestamp t0fit((AliTimestamp)(*fEvt));
436    t0fit.Add(0,0,int(t0));
437
438    // Enter the fit result as a track in the event structure
439    ntkreco++;
440    fTkfit->SetId(ntkreco);
441    fTkfit->SetCharge(fCharge);
442    fTkfit->SetParentTrack(track);
443    pos.SetTimestamp(t0fit);
444    fTkfit->SetTimestamp(t0fit);
445    fTkfit->SetReferencePoint(pos);
446    fTkfit->Set3Momentum(p);
447    for (Int_t ihit=0; ihit<fHits->GetEntries(); ihit++)
448    {
449     AliSignal* sx=(AliSignal*)fHits->At(ihit);
450     if (sx) fTkfit->AddSignal(*sx);
451    }
452    psi=GetPsi(fTkfit);
453    fFitstats->SetSignal(psi,6);
454    fFitstats->SetSignal(fPsistats.GetMedian(1),7);
455    fFitstats->SetSignal(fPsistats.GetSpread(1),8);
456    fFitstats->SetSignal(fPsistats.GetMean(1),9);
457    fFitstats->SetSignal(fPsistats.GetSigma(1),10);
458    fTkfit->SetFitDetails(fFitstats);
459    fEvt->AddTrack(fTkfit);
460   } // End loop over tracks
461  } // End loop over first guess classes
462
463 }
464 ///////////////////////////////////////////////////////////////////////////
465 void IceChi2::SetPrintLevel(Int_t level)
466 {
467 // Set the fitter (Minuit) print level.
468 // See the TFitter and TMinuit docs for details.
469 //
470 // Note : level=-2 suppresses also all fit processor warnings.
471 //        
472 // The default in the constructor is level=-2. 
473
474  fPrint=level;
475 }
476 ///////////////////////////////////////////////////////////////////////////
477 void IceChi2::UseTracks(TString classname,Int_t n)
478 {
479 // Specification of the first guess tracks to be used.
480 //
481 // classname : Specifies the first guess algorithm (e.g. "IceDwalk");
482 // n : Specifies the max. number of these tracks to be used
483 //
484 // Note : n<0 will use all the existing tracks of the specified classname
485 //
486 // The default is n=-1.
487 //
488 // Consecutive invokations of this memberfunction with different classnames
489 // will result in an incremental effect.
490 //
491 // Example :
492 // ---------
493 // UseTracks("IceDwalk",5);
494 // UseTracks("IceLinefit",2);
495 // UseTracks("IceJams");
496 //
497 // This will use the first 5 IceDwalk, the first 2 IceLinefit and all the
498 // IceJams tracks which are encountered in the event structure.
499
500  if (!fUseNames)
501  {
502   fUseNames=new TObjArray();
503   fUseNames->SetOwner();
504  }
505  
506  if (!fUseNtk) fUseNtk=new TArrayI();
507
508  // Check if this classname has already been specified before 
509  TString s;
510  Int_t nen=fUseNames->GetEntries();
511  for (Int_t i=0; i<nen; i++)
512  {
513   TObjString* sx=(TObjString*)fUseNames->At(i);
514   if (!sx) continue;
515   s=sx->GetString();
516   if (s==classname) return;
517  }
518
519  // New classname to be added into the storage
520  if (nen >= fUseNames->GetSize()) fUseNames->Expand(nen+1);
521  if (nen >= fUseNtk->GetSize()) fUseNtk->Set(nen+1);
522  
523  TObjString* name=new TObjString();
524  name->SetString(classname);
525  fUseNames->Add(name);
526  fUseNtk->AddAt(n,nen);
527 }
528 ///////////////////////////////////////////////////////////////////////////
529 void IceChi2::SelectHits(Int_t mode)
530 {
531 // Specification of the hits to be used in the minimisation.
532 //
533 // mode = 0 : All hit cleaning survived hits of the complete event are used
534 //        1 : Only the associated hits are used for each first guess track
535 //
536 // The default is mode=1.
537
538  if (mode==0 || mode==1) fSelhits=mode;
539 }
540 ///////////////////////////////////////////////////////////////////////////
541 void IceChi2::SetVgroupUsage(Int_t flag)
542 {
543 // (De)activate the distinction between v_phase and v_group of the Cherenkov light.
544 //
545 // flag = 0 : No distinction between v_phase and v_group
546 //      = 1 : Separate treatment of v_phase and v_group
547 //
548 // By default the distinction between v_phase and v_group is activated
549 // in the constructor of this class.
550  fVgroup=flag;
551 }
552 ///////////////////////////////////////////////////////////////////////////
553 void IceChi2::SetTrackName(TString s)
554 {
555 // Set (alternative) name identifier for the produced tracks.
556 // This allows unique identification of (newly) produced pandel tracks
557 // in case of re-processing of existing data with different criteria.
558 // By default the produced tracks have the name "IceChi2" which is
559 // set in the constructor of this class.
560  fTrackname=s;
561 }
562 ///////////////////////////////////////////////////////////////////////////
563 void IceChi2::SetCharge(Float_t charge)
564 {
565 // Set user defined charge for the produced tracks.
566 // This allows identification of these tracks on color displays.
567 // By default the produced tracks have charge=0 which is set in the
568 // constructor of this class.
569  fCharge=charge;
570 }
571 ///////////////////////////////////////////////////////////////////////////
572 void IceChi2::SetPenalty(Float_t val)
573 {
574 // Set user defined psi penalty value (in dB) for distance-time points that
575 // fall outside the validity rectangle.
576 // This allows investigation/tuning of the sensitivity to hits with
577 // extreme distance and/or time residual values.
578 // By default the penalty val=0 is set in the constructor of this class.
579  fPenalty=val;
580 }
581 ///////////////////////////////////////////////////////////////////////////
582 void IceChi2::FitFCN(Int_t&,Double_t*,Double_t& f,Double_t* x,Int_t)
583 {
584 // The chi-squared function used for the minimisation process.
585
586  const Float_t c=0.299792458;        // Light speed in vacuum in meters per ns
587  const Float_t npice=1.31768387;     // Phase refractive index (c/v_phase) of ice
588  const Float_t ngice=1.35075806;     // Group refractive index (c/v_group) of ice
589  const Float_t thetac=acos(1./npice);// Cherenkov angle (in radians)
590  const Double_t pi=acos(-1.);
591
592  // Angular reduction of complement of thetac due to v_phase and v_group difference
593  Float_t alphac=0;
594  if (fVgroup) alphac=atan((1.-npice/ngice)/sqrt(npice*npice-1.));
595
596  // Assumed PMT timing jitter in ns
597  const Double_t sigt=10;
598
599  f=0;
600
601  // The new r0 and p vectors and t0 from the minimisation
602  Float_t vec[3];
603  
604  AliPosition r0;
605  vec[0]=x[0];
606  vec[1]=x[1];
607  vec[2]=x[2];
608  r0.SetPosition(vec,"car");
609
610  Ali3Vector p;
611  vec[0]=1;
612  vec[1]=x[3];
613  vec[2]=x[4];
614  p.SetVector(vec,"sph");
615
616  Float_t t0=x[5];
617
618  // Construct a track with the new values from the minimisation
619  fTkfit->SetReferencePoint(r0);
620  fTkfit->Set3Momentum(p);
621
622  Int_t nhits=fHits->GetEntries();
623  AliPosition rhit;
624  Ali3Vector r12;
625  Float_t d,dist,thit,tgeo;
626  Double_t tres,chi2;
627  for (Int_t i=0; i<nhits; i++)
628  {
629   AliSignal* sx=(AliSignal*)fHits->At(i);
630   if (!sx) continue;
631   IceGOM* omx=(IceGOM*)sx->GetDevice();
632   if (!omx) continue;
633   rhit=omx->GetPosition();
634   d=fTkfit->GetDistance(rhit);
635   r12=rhit-r0;
636   dist=p.Dot(r12)+d/tan(pi/2.-thetac-alphac);
637   tgeo=t0+dist/c;
638   thit=sx->GetSignal("LE",7);
639   tres=thit-tgeo;
640
641   // Chisquared contribution for this hit
642   chi2=pow(tres/sigt,2);
643
644   f+=chi2;
645  }
646 }
647 ///////////////////////////////////////////////////////////////////////////
648 Double_t IceChi2::GetPsi(AliTrack* t)
649 {
650 // Provide Bayesian psi value for a track w.r.t. a Convoluted Pandel PDF.
651 // The Baysian psi value is defined as -loglikelihood in a decibel scale.
652 // This implies psi=-10*log10(L) where L=p(D|HI) being the likelihood of
653 // the data D under the hypothesis H and prior information I.
654 // In case of error or incomplete information a psi value of -1 is returned.
655
656  const Float_t c=0.299792458;        // Light speed in vacuum in meters per ns
657  const Float_t npice=1.31768387;     // Phase refractive index (c/v_phase) of ice
658  const Float_t ngice=1.35075806;     // Group refractive index (c/v_group) of ice
659  const Float_t thetac=acos(1./npice);// Cherenkov angle (in radians)
660  const Float_t lambda=33.3;          // Light scattering length in ice
661  const Float_t labs=98;              // Light absorbtion length in ice
662  const Float_t cice=c/ngice;         // Light speed in ice in meters per ns
663  const Float_t tau=557;
664  const Double_t rho=((1./tau)+(cice/labs));
665  const Double_t pi=acos(-1.);
666
667  // Angular reduction of complement of thetac due to v_phase and v_group difference
668  Float_t alphac=0;
669  if (fVgroup) alphac=atan((1.-npice/ngice)/sqrt(npice*npice-1.));
670
671  // Assumed PMT timing jitter in ns
672  const Double_t sigma=10;
673
674  Double_t psi=-1;
675
676  if (!t) return psi;
677
678  // The r0 and p vectors from the track
679  AliPosition* ref=t->GetReferencePoint();
680  Ali3Vector p=t->Get3Momentum();
681
682  if (!ref || p.GetNorm()<=0) return psi;
683
684  // The number of associated hits and t0 of the track
685  Int_t nhits=t->GetNsignals();
686  AliTimestamp* tstamp=ref->GetTimestamp();
687
688  if (!nhits || !tstamp) return psi;
689
690  AliPosition r0=ref->GetPosition();
691  Float_t t0=fEvt->GetDifference(tstamp,"ns");
692
693  AliPosition rhit;
694  Ali3Vector r12;
695  Float_t d,dist,thit,tgeo;
696  Double_t tres,ksi,eta,pandel;
697  Double_t cpandel1,cpandel2,cpandel3,cpandel;
698  Double_t z,k,alpha,beta,phi,n1,n2,n3,u; // Function parameters for regions 3 and 4
699  psi=0;
700  Int_t ier;
701  Double_t psihit=0;
702  fPsistats.Reset();
703  for (Int_t i=1; i<=nhits; i++)
704  {
705   AliSignal* sx=t->GetSignal(i);
706   if (!sx) continue;
707   IceGOM* omx=(IceGOM*)sx->GetDevice();
708   if (!omx) continue;
709   rhit=omx->GetPosition();
710   d=t->GetDistance(rhit);
711   ksi=d/lambda;
712   r12=rhit-r0;
713   dist=p.Dot(r12)+d/tan(pi/2.-thetac-alphac);
714   tgeo=t0+dist/c;
715   thit=sx->GetSignal("LE",7);
716   tres=thit-tgeo;
717
718   // The Convoluted Pandel function evaluation
719   // For definitions of functions and validity regions, see the
720   // CPandel writeup of O. Fadiran, G. Japaridze and N. van Eijndhoven
721
722   // Move points which are outside the validity rectangle in the (tres,ksi) space
723   // to the edge of the validity rectangle and signal the use of the penalty
724   ier=0;
725   if (tres<-25.*sigma)
726   {
727    tres=-25.*sigma;
728    ier=1;
729   }
730   if (tres>3500)
731   {
732    tres=3500;
733    ier=1;
734   }
735   if (ksi>50)
736   {
737    ksi=50;
738    ier=1;
739   }
740
741   eta=(rho*sigma)-(tres/sigma);
742
743   if (ksi<=0) // The zero distance (ksi=0) axis
744   {
745    cpandel=exp(-tres*tres/(2.*sigma*sigma))/(sigma*sqrt(2.*pi));
746   }
747   else if (ksi<=5 && tres>=-5.*sigma && tres<=30.*sigma) // The exact expression in region 1
748   {
749    cpandel1=pow(rho,ksi)*pow(sigma,ksi-1.)*exp(-tres*tres/(2.*sigma*sigma))/pow(2.,0.5*(1.+ksi));
750    cpandel2=ROOT::Math::conf_hyperg(ksi/2.,0.5,eta*eta/2.)/TMath::Gamma((ksi+1.)/2.);
751    cpandel3=sqrt(2.)*eta*ROOT::Math::conf_hyperg((ksi+1.)/2.,1.5,eta*eta/2.)/TMath::Gamma(ksi/2.);
752
753    cpandel=cpandel1*(cpandel2-cpandel3);
754   }
755   else if (ksi<=1 && tres>30.*sigma && tres<=3500) // Approximation in region 2
756   {
757    pandel=pow(rho,ksi)*pow(tres,(ksi-1.))*exp(-rho*tres)/TMath::Gamma(ksi);
758
759    cpandel=exp(rho*rho*sigma*sigma/2.)*pandel;
760   }
761   else if (ksi<=1 && tres<-5.*sigma && tres>=-25.*sigma) // Approximation in region 5
762   {
763    cpandel=pow(rho*sigma,ksi)*pow(eta,-ksi)*exp(-tres*tres/(2.*sigma*sigma))/(sigma*sqrt(2.*pi));
764   }
765   else if (ksi<=50 && tres>=0 && tres<=3500) // Approximation in region 3
766   {
767    z=-eta/sqrt(4.*ksi-2.);
768    k=0.5*(z*sqrt(1.+z*z)+log(z+sqrt(1.+z*z)));
769    alpha=-tres*tres/(2.*sigma*sigma)+eta*eta/4.-ksi/2.+0.25+k*(2.*ksi-1.);
770    alpha+=-log(1.+z*z)/4.-ksi*log(2.)/2.+(ksi-1.)*log(2.*ksi-1.)/2.+ksi*log(rho)+(ksi-1.)*log(sigma);
771    beta=0.5*(z/sqrt(1.+z*z)-1.);
772    n1=beta*(20.*beta*beta+30.*beta+9.)/12.;
773    n2=pow(beta,2.)*(6160.*pow(beta,4.)+18480.*pow(beta,3.)+19404.*pow(beta,2.)+8028.*beta+945.)/288.;
774    n3=27227200.*pow(beta,6.)+122522400.*pow(beta,5.)+220540320.*pow(beta,4.);
775    n3+=200166120.*pow(beta,3.)+94064328.*pow(beta,2.)+20546550.*beta+1403325.;
776    n3*=pow(beta,3.)/51840.;
777    phi=1.-n1/(2.*ksi-1.)+n2/pow(2.*ksi-1.,2.)-n3/pow(2.*ksi-1.,3.);
778    cpandel=exp(alpha)*phi/TMath::Gamma(ksi);
779   }
780   else if (ksi<=50 && tres<0 && tres>=-25.*sigma) // Approximation in region 4
781   {
782    z=eta/sqrt(4.*ksi-2.);
783    k=0.5*(z*sqrt(1.+z*z)+log(z+sqrt(1.+z*z)));
784    u=exp(ksi/2.-0.25)*pow(2.*ksi-1.,-ksi/2.)*pow(2.,(ksi-1.)/2.);
785    beta=0.5*(z/sqrt(1.+z*z)-1.);
786    n1=beta*(20.*beta*beta+30.*beta+9.)/12.;
787    n2=pow(beta,2.)*(6160.*pow(beta,4.)+18480.*pow(beta,3.)+19404.*pow(beta,2.)+8028.*beta+945.)/288.;
788    n3=27227200.*pow(beta,6.)+122522400.*pow(beta,5.)+220540320.*pow(beta,4.);
789    n3+=200166120.*pow(beta,3.)+94064328.*pow(beta,2.)+20546550.*beta+1403325.;
790    n3*=pow(beta,3.)/51840.;
791    phi=1.+n1/(2.*ksi-1.)+n2/pow(2.*ksi-1.,2.)+n3/pow(2.*ksi-1.,3.);
792    cpandel=pow(rho,ksi)*pow(sigma,ksi-1.)*exp(-pow(tres,2.)/(2.*pow(sigma,2.))+pow(eta,2.)/4.)/sqrt(2.*pi);
793    cpandel*=u*phi*exp(-k*(2.*ksi-1.))*pow(1.+z*z,-0.25);
794   }
795   else // (tres,ksi) outside validity rectangle
796   {
797    ier=1;
798    cout << " *IceChi2::GetPsi* Outside rectangle. We should never get here." << endl;
799   }
800
801   // Use 10*log10 expression to obtain intuitive dB scale
802   // Omit (small) negative values which are possible due to computer accuracy
803   psihit=0;
804   if (cpandel>0) psihit=-10.*log10(cpandel);
805
806   // Penalty in dB for (tres,ksi) points outside the validity rectangle
807   if (ier) psihit+=fPenalty;
808
809   // Update the psi statistics for this hit
810   fPsistats.Enter(float(psihit));
811   psi+=psihit;
812  }
813  return psi;
814 }
815 ///////////////////////////////////////////////////////////////////////////