1 /*******************************************************************************
2 * Copyright(c) 2003, IceCube Experiment at the South Pole. All rights reserved.
4 * Author: The IceCube RALICE-based Offline Project.
5 * Contributors are mentioned in the code where appropriate.
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 *******************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////
20 // TTask derived class to perform track fitting via minimisation of a
21 // convoluted Pandel pdf.
23 // The code in this processor is based on the algorithms as developed
24 // by Oladipo Fadiran, George Japaridze (Clark Atlanta University, USA)
25 // and Nick van Eijndhoven (Utrecht University, The Netherlands).
27 // For the minimisation process the TFitter facility, which is basically Minuit,
28 // is used. Minimisation is performed by invokation of the SIMPLEX method,
29 // followed by an invokation of HESSE to determine the uncertainties on the results.
30 // The statistics of the TFitter result are stored as an AliSignal object
31 // in the track, which can be obtained via the GetFitDetails memberfunction.
32 // In the minimisation procedure an overall plausibility for the fitted track
33 // is determined based on a convoluted Pandel pdf value for each used hit.
34 // This track plausibility is expressed in terms of a Bayesian psi value
35 // w.r.t. a Convoluted Pandel PDF.
36 // The Baysian psi value is defined as -loglikelihood in a decibel scale.
37 // This implies psi=-10*log10(L) where L=p(D|HI) being the likelihood of
38 // the data D under the hypothesis H and prior information I.
39 // Since all (associated) hits contribute independently to the Bayesian psi
40 // value, this psi value is built up by summation of the various hit contributions.
41 // As such, the FitDetails entries contain the statistics of all the different
42 // hit contributions, like PsiMedian, PsiMean, and PsiSigma.
43 // The Bayesian psi value is available in the fit details under the name "PsiSum".
44 // In addition the standard Minuit results like IERFIT, FCN, EDM etc... are
45 // also available from the FitDetails.
47 // The convoluted Pandel value is evaluated in various areas in the distance-time
48 // space as described in the CPandel writeup by O. Fadiran, G. Japaridze
49 // and N. van Eijndhoven.
50 // In case the distance-time point of a certain hit falls outside the
51 // validity rectangle, the point is moved onto the corresponding side location
52 // of the rectangle. For this new location the Pandel value is evaluated for
53 // this hit and an extra penalty is added to the corresponding psi value
55 // By default this penalty value amounts to 0 dB, but the user can
56 // modify this penalty value via the memberfunction SetPenalty.
57 // This allows investigation/tuning of the sensitivity to hits with
58 // extreme distance and/or time residual values.
60 // A separate treatment of the phase and group velocities is introduced
61 // which will provide more accurate time residuals due to the different
62 // velocities of the Cerenkov wave front (v_phase) and the actually detected
64 // This distinction between v_phase and v_group can be (de)activated via the
65 // memberfunction SetVgroupUsage(). By default the distinction between v_phase
66 // and v_group is activated in the constructor of this class.
68 // Use the UseTracks memberfunction to specify the first guess tracks
69 // to be processed by the minimiser.
70 // By default only the first encountered IceDwalk track will be processed.
72 // Use the SelectHits memberfunction to specify the hits to be used.
73 // By default only the hits associated to the first guess track are used.
75 // Information about the actual parameter settings can be found in the event
76 // structure itself via the device named "IcePandel".
78 // The fit processor printlevel can be selected via the memberfunction SetPrintLevel.
79 // By default all printout is suppressed (i.e. level=-2).
81 // An example of how to invoke this processor after Xtalk, hit cleaning
82 // and a direct walk first guess estimate can be found in the ROOT example
83 // macro icepandel.cc which resides in the /macros subdirectory.
85 // The minimisation results are stored in the IceEvent structure as
86 // tracks with as default the name "IcePandel" (just like the first guess
87 // results of e.g. IceDwalk).
88 // This track name identifier can be modified by the user via the
89 // SetTrackName() memberfunction. This will allow unique identification
90 // of tracks which are produced when re-processing existing data with
91 // different criteria.
92 // By default the charge of the produced tracks is set to 0, since
93 // no distinction can be made between positive or negative tracks.
94 // However, the user can define the track charge by invokation
95 // of the memberfunction SetCharge().
96 // This facility may be used to distinguish tracks produced by the
97 // various reconstruction algorithms in a (3D) colour display
98 // (see the class AliHelix for further details).
99 // A pointer to the first guess track which was used as input is available
100 // via the GetParentTrack facility of these "IcePandel" tracks.
101 // Furthermore, all the hits that were used in the minisation are available
102 // via the GetSignal facility of a certain track.
104 // An example of how the various data can be accessed is given below,
105 // where "evt" indicates the pointer to the IceEvent structure.
107 // Example for accessing data :
108 // ----------------------------
109 // TObjArray* tracks=evt->GetTracks("IcePandel");
110 // if (!tracks) return;
111 // AliPosition* r0=0;
113 // for (Int_t jtk=0; jtk<tracks->GetEntries(); jtk++)
115 // AliTrack* tx=(AliTrack*)tracks->At(jtk);
116 // if (!tx) continue;
118 // r0=tx->GetReferencePoint();
119 // if (r0) r0->Data();
120 // sx=(AliSignal*)tx->GetFitDetails();
121 // if (sx) fcn=sx->GetSignal("FCN");
122 // AliTrack* tx2=tx->GetParentTrack();
123 // if (!tx2) continue;
125 // r0=tx2->GetReferencePoint();
126 // if (r0) r0->Data();
131 // 1) This processor only works properly on data which are Time and ADC
132 // calibrated and contain tracks from first guess algorithms like
134 // 2) In view of the usage of TFitter/Minuit minimisation, a global pointer
135 // to the instance of this class (gIcePandel) and a global static
136 // wrapper function (IcePandelFCN) have been introduced, to allow the
137 // actual minimisation to be performed via the memberfunction FitFCN.
138 // This implies that in a certain processing job only 1 instance of
139 // this IcePandel class may occur.
141 //--- Author: Nick van Eijndhoven 09-feb-2006 Utrecht University
142 //- Modified: NvE $Date$ Utrecht University
143 ///////////////////////////////////////////////////////////////////////////
145 #include "IcePandel.h"
146 #include "Riostream.h"
148 // Global pointer to the instance of this object
149 IcePandel* gIcePandel=0;
151 // TFitter/Minuit interface to IcePandel::FitFCN
152 void IcePandelFCN(Int_t& npar,Double_t* gin,Double_t& f,Double_t* u,Int_t flag)
154 if (gIcePandel) gIcePandel->FitFCN(npar,gin,f,u,flag);
157 ClassImp(IcePandel) // Class implementation to enable ROOT I/O
159 IcePandel::IcePandel(const char* name,const char* title) : TTask(name,title)
161 // Default constructor.
171 fTrackname="IcePandel";
177 // Set the global pointer to this instance
180 ///////////////////////////////////////////////////////////////////////////
181 IcePandel::~IcePandel()
183 // Default destructor.
215 ///////////////////////////////////////////////////////////////////////////
216 void IcePandel::Exec(Option_t* opt)
218 // Implementation of the hit fitting procedure.
221 AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
225 fEvt=(IceEvent*)parent->GetObject("IceEvent");
228 // Storage of the used parameters in the IcePandel device
230 params.SetNameTitle("IcePandel","IcePandel processor parameters");
231 params.SetSlotName("Selhits",1);
232 params.SetSlotName("Penalty",2);
233 params.SetSlotName("Vgroup",3);
235 params.SetSignal(fSelhits,1);
236 params.SetSignal(fPenalty,2);
237 params.SetSignal(fVgroup,3);
239 fEvt->AddDevice(params);
241 if (!fUseNames) UseTracks("IceDwalk",1);
243 Int_t nclasses=fUseNames->GetEntries(); // Number of first guess classes to be processed
244 Int_t ntkmax=0; // Max. number of tracks for a certain class
250 cout << " *IcePandel* First guess selections to be processed (-1=all)." << endl;
251 for (Int_t i=0; i<nclasses; i++)
253 strx=(TObjString*)fUseNames->At(i);
255 str=strx->GetString();
256 ntkmax=fUseNtk->At(i);
257 cout << " Maximally " << ntkmax << " track(s) per event for procedure : " << str.Data() << endl;
259 cout << " *IcePandel* Hit selection mode : " << fSelhits << endl;
260 cout << " *IcePandel* Penalty value for minimiser : " << fPenalty << " dB." << endl;
263 fPsistats.SetStoreMode(1);
268 const Double_t pi=acos(-1.);
269 const Double_t e=exp(1.);
271 // Initialisation of the minimisation processor
272 Double_t arglist[100];
273 if (!fFitter) fFitter=new TFitter();
275 // The number of reconstructed tracks already present in the event
276 Int_t ntkreco=fEvt->GetNtracks(1);
280 fHits=new TObjArray();
287 // If selected, use all the good quality hits of the complete event
290 TObjArray* hits=fEvt->GetHits("IceGOM");
291 for (Int_t ih=0; ih<hits->GetEntries(); ih++)
293 AliSignal* sx=(AliSignal*)hits->At(ih);
295 if (sx->GetDeadValue("ADC") || sx->GetDeadValue("LE") || sx->GetDeadValue("TOT")) continue;
300 // Track by track processing of the selected first guess classes
307 fTkfit=new AliTrack();
308 fTkfit->SetNameTitle(fTrackname.Data(),"IcePandel fit result");
312 fFitstats=new AliSignal();
313 fFitstats->SetNameTitle("Fitstats","TFitter stats for Pandel fit");
314 fFitstats->SetSlotName("IERFIT",1);
315 fFitstats->SetSlotName("FCN",2);
316 fFitstats->SetSlotName("EDM",3);
317 fFitstats->SetSlotName("NVARS",4);
318 fFitstats->SetSlotName("IERERR",5);
319 fFitstats->SetSlotName("PsiSum",6);
320 fFitstats->SetSlotName("PsiMedian",7);
321 fFitstats->SetSlotName("PsiSpread",8);
322 fFitstats->SetSlotName("PsiMean",9);
323 fFitstats->SetSlotName("PsiSigma",10);
325 Float_t x,y,z,theta,phi,t0;
326 Double_t amin,edm,errdef; // Minimisation stats
327 Int_t ierfit,iererr,nvpar,nparx; // Minimisation stats
331 for (Int_t iclass=0; iclass<nclasses; iclass++) // Loop over first guess classes
333 strx=(TObjString*)fUseNames->At(iclass);
335 str=strx->GetString();
336 ntkmax=fUseNtk->At(iclass);
337 TObjArray* tracks=fEvt->GetTracks(str);
338 ntk=tracks->GetEntries();
339 if (ntkmax>0 && ntk>ntkmax) ntk=ntkmax;
341 for (Int_t jtk=0; jtk<ntk; jtk++) // Loop over tracks of a certain class
343 track=(AliTrack*)tracks->At(jtk);
344 if (!track) continue;
346 AliPosition* r0=track->GetReferencePoint();
349 AliTimestamp* tt0=r0->GetTimestamp();
351 // If selected, use only the first guess track associated hits
355 nsig=track->GetNsignals();
356 for (Int_t is=1; is<=nsig; is++)
358 AliSignal* sx=track->GetSignal(is);
360 if (!sx->GetDevice()->InheritsFrom("IceGOM")) continue;
361 if (sx->GetDeadValue("ADC") || sx->GetDeadValue("LE") || sx->GetDeadValue("TOT")) continue;
366 if (!fHits->GetEntries()) continue;
368 r0->GetVector(vec,"car");
369 r0->GetErrors(err,"car");
375 p=track->Get3Momentum();
376 p.GetVector(vec,"sph");
381 t0=fEvt->GetDifference(tt0,"ns");
383 // Process this first guess track with its associated hits
386 // Set user selected TFitter printout level
388 if (fPrint==-2) arglist[0]=-1;
389 fFitter->ExecuteCommand("SET PRINT",arglist,1);
390 if (fPrint==-2) fFitter->ExecuteCommand("SET NOWARNINGS",arglist,0);
392 fFitter->SetFitMethod("loglikelihood");
394 // Define errors to represent 1 sigma for this likelihood scale
395 arglist[0]=5.*log10(e);
396 fFitter->ExecuteCommand("SET ERRORDEF",arglist,1);
398 fFitter->SetParameter(0,"r0x",x,0.1,0,0);
399 fFitter->SetParameter(1,"r0y",y,0.1,0,0);
400 fFitter->SetParameter(2,"r0z",z,0.1,0,0);
401 fFitter->SetParameter(3,"theta",theta,0.001,0,pi);
402 fFitter->SetParameter(4,"phi",phi,0.001,0,2.*pi);
403 fFitter->SetParameter(5,"t0",t0,1.,0,32000);
405 fFitter->SetFCN(IcePandelFCN);
410 ierfit=fFitter->ExecuteCommand("SIMPLEX",arglist,0);
412 fFitter->GetStats(amin,edm,errdef,nvpar,nparx);
415 fFitstats->SetSignal(ierfit,1);
416 fFitstats->SetSignal(amin,2);
417 fFitstats->SetSignal(edm,3);
418 fFitstats->SetSignal(nvpar,4);
420 fFitstats->SetSignal(fPsistats.GetSum(1),6);
421 fFitstats->SetSignal(fPsistats.GetMedian(1),7);
422 fFitstats->SetSignal(fPsistats.GetSpread(1),8);
423 fFitstats->SetSignal(fPsistats.GetMean(1),9);
424 fFitstats->SetSignal(fPsistats.GetSigma(1),10);
426 iererr=fFitter->ExecuteCommand("HESSE",arglist,0);
427 fFitstats->SetSignal(iererr,5);
429 // Resulting parameters after minimisation and error calculation
430 vec[0]=fFitter->GetParameter(0);
431 vec[1]=fFitter->GetParameter(1);
432 vec[2]=fFitter->GetParameter(2);
433 err[0]=fFitter->GetParError(0);
434 err[1]=fFitter->GetParError(1);
435 err[2]=fFitter->GetParError(2);
436 pos.SetPosition(vec,"car");
437 pos.SetPositionErrors(err,"car");
440 vec[1]=fFitter->GetParameter(3);
441 vec[2]=fFitter->GetParameter(4);
443 err[1]=fFitter->GetParError(3);
444 err[2]=fFitter->GetParError(4);
445 p.SetVector(vec,"sph");
446 p.SetErrors(err,"sph");
448 t0=fFitter->GetParameter(5);
449 AliTimestamp t0fit((AliTimestamp)(*fEvt));
450 t0fit.Add(0,0,int(t0));
452 // Enter the fit result as a track in the event structure
454 fTkfit->SetId(ntkreco);
455 fTkfit->SetCharge(fCharge);
456 fTkfit->SetParentTrack(track);
457 pos.SetTimestamp(t0fit);
458 fTkfit->SetTimestamp(t0fit);
459 fTkfit->SetReferencePoint(pos);
460 fTkfit->Set3Momentum(p);
461 for (Int_t ihit=0; ihit<fHits->GetEntries(); ihit++)
463 AliSignal* sx=(AliSignal*)fHits->At(ihit);
464 if (sx) fTkfit->AddSignal(*sx);
466 fTkfit->SetFitDetails(fFitstats);
467 fEvt->AddTrack(fTkfit);
468 } // End loop over tracks
469 } // End loop over first guess classes
472 ///////////////////////////////////////////////////////////////////////////
473 void IcePandel::SetPrintLevel(Int_t level)
475 // Set the fitter (Minuit) print level.
476 // See the TFitter and TMinuit docs for details.
478 // Note : level=-2 suppresses also all fit processor warnings.
480 // The default in the constructor is level=-2.
484 ///////////////////////////////////////////////////////////////////////////
485 void IcePandel::UseTracks(TString classname,Int_t n)
487 // Specification of the first guess tracks to be used.
489 // classname : Specifies the first guess algorithm (e.g. "IceDwalk");
490 // n : Specifies the max. number of these tracks to be used
492 // Note : n<0 will use all the existing tracks of the specified classname
494 // The default is n=-1.
496 // Consecutive invokations of this memberfunction with different classnames
497 // will result in an incremental effect.
501 // UseTracks("IceDwalk",5);
502 // UseTracks("IceLinefit",2);
503 // UseTracks("IceJams");
505 // This will use the first 5 IceDwalk, the first 2 IceLinefit and all the
506 // IceJams tracks which are encountered in the event structure.
510 fUseNames=new TObjArray();
511 fUseNames->SetOwner();
514 if (!fUseNtk) fUseNtk=new TArrayI();
516 // Check if this classname has already been specified before
518 Int_t nen=fUseNames->GetEntries();
519 for (Int_t i=0; i<nen; i++)
521 TObjString* sx=(TObjString*)fUseNames->At(i);
524 if (s==classname) return;
527 // New classname to be added into the storage
528 if (nen >= fUseNames->GetSize()) fUseNames->Expand(nen+1);
529 if (nen >= fUseNtk->GetSize()) fUseNtk->Set(nen+1);
531 TObjString* name=new TObjString();
532 name->SetString(classname);
533 fUseNames->Add(name);
534 fUseNtk->AddAt(n,nen);
536 ///////////////////////////////////////////////////////////////////////////
537 void IcePandel::SelectHits(Int_t mode)
539 // Specification of the hits to be used in the minimisation.
541 // mode = 0 : All hit cleaning survived hits of the complete event are used
542 // 1 : Only the associated hits are used for each first guess track
544 // The default is mode=1.
546 if (mode==0 || mode==1) fSelhits=mode;
548 ///////////////////////////////////////////////////////////////////////////
549 void IcePandel::SetVgroupUsage(Int_t flag)
551 // (De)activate the distinction between v_phase and v_group of the Cherenkov light.
553 // flag = 0 : No distinction between v_phase and v_group
554 // = 1 : Separate treatment of v_phase and v_group
556 // By default the distinction between v_phase and v_group is activated
557 // in the constructor of this class.
560 ///////////////////////////////////////////////////////////////////////////
561 void IcePandel::SetTrackName(TString s)
563 // Set (alternative) name identifier for the produced tracks.
564 // This allows unique identification of (newly) produced pandel tracks
565 // in case of re-processing of existing data with different criteria.
566 // By default the produced tracks have the name "IcePandel" which is
567 // set in the constructor of this class.
570 ///////////////////////////////////////////////////////////////////////////
571 void IcePandel::SetCharge(Float_t charge)
573 // Set user defined charge for the produced tracks.
574 // This allows identification of these tracks on color displays.
575 // By default the produced tracks have charge=0 which is set in the
576 // constructor of this class.
579 ///////////////////////////////////////////////////////////////////////////
580 void IcePandel::SetPenalty(Float_t val)
582 // Set user defined psi penalty value (in dB) in the minimiser for
583 // distance-time points that fall outside the validity rectangle.
584 // This allows investigation/tuning of the sensitivity to hits with
585 // extreme distance and/or time residual values.
586 // By default the penalty val=0 is set in the constructor of this class.
589 ///////////////////////////////////////////////////////////////////////////
590 void IcePandel::FitFCN(Int_t&,Double_t*,Double_t& f,Double_t* x,Int_t)
592 // Minimisation of the Bayesian psi value for a track w.r.t. a Convoluted Pandel PDF.
593 // The Baysian psi value is defined as -loglikelihood in a decibel scale.
594 // This implies psi=-10*log10(L) where L=p(D|HI) being the likelihood of
595 // the data D under the hypothesis H and prior information I.
597 const Float_t c=0.299792458; // Light speed in vacuum in meters per ns
598 const Float_t npice=1.31768387; // Phase refractive index (c/v_phase) of ice
599 const Float_t ngice=1.35075806; // Group refractive index (c/v_group) of ice
600 const Float_t thetac=acos(1./npice);// Cherenkov angle (in radians)
601 const Float_t lambda=33.3; // Light scattering length in ice
602 const Float_t labs=98; // Light absorbtion length in ice
603 const Float_t cice=c/ngice; // Light speed in ice in meters per ns
604 const Float_t tau=557;
605 const Double_t rho=((1./tau)+(cice/labs));
606 const Double_t pi=acos(-1.);
608 // Angular reduction of complement of thetac due to v_phase and v_group difference
610 if (fVgroup) alphac=atan((1.-npice/ngice)/sqrt(npice*npice-1.));
612 // Assumed PMT timing jitter in ns
613 const Double_t sigma=10;
617 // The new r0 and p vectors and t0 from the minimisation
624 r0.SetPosition(vec,"car");
630 p.SetVector(vec,"sph");
634 // Construct a track with the new values from the minimisation
635 fTkfit->SetReferencePoint(r0);
636 fTkfit->Set3Momentum(p);
638 Int_t nhits=fHits->GetEntries();
642 Float_t d,dist,thit,tgeo;
643 Double_t tres,ksi,eta,pandel;
644 Double_t cpandel1,cpandel2,cpandel3,cpandel;
645 Double_t z,k,alpha,beta,phi,n1,n2,n3,u; // Function parameters for regions 3 and 4
649 //@@@ for (Int_t i=1; i<=nhits; i++)
650 for (Int_t i=0; i<nhits; i++)
652 AliSignal* sx=(AliSignal*)fHits->At(i);
654 IceGOM* omx=(IceGOM*)sx->GetDevice();
656 rhit=omx->GetPosition();
657 d=fTkfit->GetDistance(rhit);
660 dist=p.Dot(r12)+d/tan(pi/2.-thetac-alphac);
662 thit=sx->GetSignal("LE",7);
665 // The Convoluted Pandel function evaluation
666 // For definitions of functions and validity regions, see the
667 // CPandel writeup of O. Fadiran, G. Japaridze and N. van Eijndhoven
669 // Move points which are outside the validity rectangle in the (tres,ksi) space
670 // to the edge of the validity rectangle and signal the use of the penalty
688 eta=(rho*sigma)-(tres/sigma);
690 if (ksi<=0) // The zero distance (ksi=0) axis
692 cpandel=exp(-tres*tres/(2.*sigma*sigma))/(sigma*sqrt(2.*pi));
694 else if (ksi<=5 && tres>=-5.*sigma && tres<=30.*sigma) // The exact expression in region 1
696 cpandel1=pow(rho,ksi)*pow(sigma,ksi-1.)*exp(-tres*tres/(2.*sigma*sigma))/pow(2.,0.5*(1.+ksi));
697 cpandel2=ROOT::Math::conf_hyperg(ksi/2.,0.5,eta*eta/2.)/TMath::Gamma((ksi+1.)/2.);
698 cpandel3=sqrt(2.)*eta*ROOT::Math::conf_hyperg((ksi+1.)/2.,1.5,eta*eta/2.)/TMath::Gamma(ksi/2.);
700 cpandel=cpandel1*(cpandel2-cpandel3);
702 else if (ksi<=1 && tres>30.*sigma && tres<=3500) // Approximation in region 2
704 pandel=pow(rho,ksi)*pow(tres,(ksi-1.))*exp(-rho*tres)/TMath::Gamma(ksi);
706 cpandel=exp(rho*rho*sigma*sigma/2.)*pandel;
708 else if (ksi<=1 && tres<-5.*sigma && tres>=-25.*sigma) // Approximation in region 5
710 cpandel=pow(rho*sigma,ksi)*pow(eta,-ksi)*exp(-tres*tres/(2.*sigma*sigma))/(sigma*sqrt(2.*pi));
712 else if (ksi<=50 && tres>=0 && tres<=3500) // Approximation in region 3
714 z=-eta/sqrt(4.*ksi-2.);
715 k=0.5*(z*sqrt(1.+z*z)+log(z+sqrt(1.+z*z)));
716 alpha=-tres*tres/(2.*sigma*sigma)+eta*eta/4.-ksi/2.+0.25+k*(2.*ksi-1.);
717 alpha+=-log(1.+z*z)/4.-ksi*log(2.)/2.+(ksi-1.)*log(2.*ksi-1.)/2.+ksi*log(rho)+(ksi-1.)*log(sigma);
718 beta=0.5*(z/sqrt(1.+z*z)-1.);
719 n1=beta*(20.*beta*beta+30.*beta+9.)/12.;
720 n2=pow(beta,2.)*(6160.*pow(beta,4.)+18480.*pow(beta,3.)+19404.*pow(beta,2.)+8028.*beta+945.)/288.;
721 n3=27227200.*pow(beta,6.)+122522400.*pow(beta,5.)+220540320.*pow(beta,4.);
722 n3+=200166120.*pow(beta,3.)+94064328.*pow(beta,2.)+20546550.*beta+1403325.;
723 n3*=pow(beta,3.)/51840.;
724 phi=1.-n1/(2.*ksi-1.)+n2/pow(2.*ksi-1.,2.)-n3/pow(2.*ksi-1.,3.);
725 cpandel=exp(alpha)*phi/TMath::Gamma(ksi);
727 else if (ksi<=50 && tres<0 && tres>=-25.*sigma) // Approximation in region 4
729 z=eta/sqrt(4.*ksi-2.);
730 k=0.5*(z*sqrt(1.+z*z)+log(z+sqrt(1.+z*z)));
731 u=exp(ksi/2.-0.25)*pow(2.*ksi-1.,-ksi/2.)*pow(2.,(ksi-1.)/2.);
732 beta=0.5*(z/sqrt(1.+z*z)-1.);
733 n1=beta*(20.*beta*beta+30.*beta+9.)/12.;
734 n2=pow(beta,2.)*(6160.*pow(beta,4.)+18480.*pow(beta,3.)+19404.*pow(beta,2.)+8028.*beta+945.)/288.;
735 n3=27227200.*pow(beta,6.)+122522400.*pow(beta,5.)+220540320.*pow(beta,4.);
736 n3+=200166120.*pow(beta,3.)+94064328.*pow(beta,2.)+20546550.*beta+1403325.;
737 n3*=pow(beta,3.)/51840.;
738 phi=1.+n1/(2.*ksi-1.)+n2/pow(2.*ksi-1.,2.)+n3/pow(2.*ksi-1.,3.);
739 cpandel=pow(rho,ksi)*pow(sigma,ksi-1.)*exp(-pow(tres,2.)/(2.*pow(sigma,2.))+pow(eta,2.)/4.)/sqrt(2.*pi);
740 cpandel*=u*phi*exp(-k*(2.*ksi-1.))*pow(1.+z*z,-0.25);
742 else // (tres,ksi) outside validity rectangle
745 cout << " *IcePandel::FitFCN* Outside rectangle. We should never get here." << endl;
748 // Use 10*log10 expression to obtain intuitive dB scale
749 // Omit (small) negative values which are possible due to computer accuracy
751 if (cpandel>0) psihit=-10.*log10(cpandel);
753 // Penalty in dB for (tres,ksi) points outside the validity rectangle
754 if (ier) psihit+=fPenalty;
756 // Update the psi statistics for this hit
757 fPsistats.Enter(float(psihit));
761 ///////////////////////////////////////////////////////////////////////////