]>
Commit | Line | Data |
---|---|---|
219e9b7e | 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 | // | |
25eefd00 | 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 | // | |
219e9b7e | 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 | // | |
25eefd00 | 71 | // Information about the actual parameter settings can be found in the event |
72 | // structure itself via the device named "IceChi2". | |
73 | // | |
219e9b7e | 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; | |
25eefd00 | 161 | fVgroup=1; |
219e9b7e | 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 | ||
25eefd00 | 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 | ||
219e9b7e | 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); | |
25eefd00 | 316 | fFitstats->SetSlotName("PsiSpread",8); |
317 | fFitstats->SetSlotName("PsiMean",9); | |
318 | fFitstats->SetSlotName("PsiSigma",10); | |
219e9b7e | 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)); | |
b57fed05 | 436 | t0fit.Add(0,0,int(t0)); |
219e9b7e | 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); | |
25eefd00 | 455 | fFitstats->SetSignal(fPsistats.GetSpread(1),8); |
456 | fFitstats->SetSignal(fPsistats.GetMean(1),9); | |
457 | fFitstats->SetSignal(fPsistats.GetSigma(1),10); | |
219e9b7e | 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 | /////////////////////////////////////////////////////////////////////////// | |
25eefd00 | 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 | /////////////////////////////////////////////////////////////////////////// | |
219e9b7e | 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 | ||
25eefd00 | 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.)); | |
219e9b7e | 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; | |
25eefd00 | 636 | dist=p.Dot(r12)+d/tan(pi/2.-thetac-alphac); |
219e9b7e | 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 | ||
25eefd00 | 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) | |
219e9b7e | 660 | const Float_t lambda=33.3; // Light scattering length in ice |
661 | const Float_t labs=98; // Light absorbtion length in ice | |
25eefd00 | 662 | const Float_t cice=c/ngice; // Light speed in ice in meters per ns |
219e9b7e | 663 | const Float_t tau=557; |
664 | const Double_t rho=((1./tau)+(cice/labs)); | |
665 | const Double_t pi=acos(-1.); | |
666 | ||
25eefd00 | 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 | ||
219e9b7e | 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; | |
25eefd00 | 713 | dist=p.Dot(r12)+d/tan(pi/2.-thetac-alphac); |
219e9b7e | 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 | /////////////////////////////////////////////////////////////////////////// |