]>
Commit | Line | Data |
---|---|---|
25eefd00 | 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 IceDwalkx | |
20 | // TTask derived class to perform (Amanda-like) direct walk track reconstruction. | |
21 | // This class is kept to provide a procedure that can closely match the | |
22 | // performance of the original Amanda direct walk-II algorithm as implemented | |
23 | // in the Sieglinde reconstruction framework. | |
24 | // For a more sophisticated direct walk reconstruction procedure, please refer | |
25 | // to the class IceDwalk. | |
26 | // The IceDwalkx procedure is based on the method described in the Amanda | |
27 | // publication in Nuclear Instruments and Methods A524 (2004) 179-180. | |
28 | // However, the Amanda method has been extended with the intention to | |
29 | // take also multiple (muon) tracks within 1 event into account. | |
30 | // This will not only provide a means to reconstruct muon bundles and | |
31 | // multiple track events in IceCube, but will also allow to reduce the | |
32 | // background of faked upgoing muons as a result of multiple downgoing | |
33 | // muons hitting the top and bottom parts of the detector. | |
34 | // A further extension of the original Amanda method is the separate treatment | |
35 | // of the phase and group velocities as introduced in collaboration with | |
36 | // George Japaridze (Clark Atlanta University, USA) which will provide more | |
37 | // accurate time residuals due to the different velocities of the Cerenkov | |
38 | // wave front (v_phase) and the actually detected photons (v_group). | |
39 | // This distinction between v_phase and v_group can be (de)activated via the | |
40 | // memberfunction SetVgroupUsage(). By default the distinction between v_phase | |
41 | // and v_group is activated in the constructor of this class. | |
42 | // To prevent waisting CPU time in trying to reconstruct (high-energy) cascade | |
43 | // events, or to select specifically reconstruction of low multiplicity events, | |
44 | // the user may invoke the memberfunctions SetMaxModA() and SetMinModA(). | |
45 | // This allows selection of events for processing with a certain maximum and/or | |
46 | // minimum number of good Amanda OMs firing. | |
47 | // By default the minimum and maximum are set to 0 and 999, respectively, | |
48 | // in the constructor, which implies no multiplicity selection. | |
49 | // The maximum number of good hits per Amanda OM to be used for the reconstruction | |
50 | // can be specified via the memberfunction SetMaxHitsA(). | |
51 | // By default only the first good hit of each Amanda OM is used to be consistent | |
52 | // with the original Sieglinde direct walk procedure of the above NIM article. | |
53 | // Note that when all the good hits of an OM are used, this may lead to large | |
54 | // processing time in case many noise and/or afterpulse signals are not | |
55 | // recognised by the hit cleaning procedure. | |
56 | // | |
57 | // Information about the actual parameter settings can be found in the event | |
58 | // structure itself via the device named "IceDwalkx". | |
59 | // | |
60 | // The various reconstruction steps are summarised as follows : | |
61 | // | |
62 | // 1) Construction of track elements (TE's). | |
63 | // A track element is a straight line connecting two hits that | |
64 | // appeared at some minimum distance d and within some maximum | |
65 | // time difference dt. | |
66 | // The default values for d and dt are given in eq. (20) of the | |
67 | // NIM article, but can be modified by the appropriate Set functions. | |
68 | // For dt a default margin of 30 ns is used (according to eq. (20)), | |
69 | // but also this margin may be modified via the appropriate Set function. | |
70 | // The reference point r0 of the TE is taken as the center between | |
71 | // the two hit positions and the TE timestamp t0 at the position r0 | |
72 | // is taken as the IceEvent timestamp increased by the average of the | |
73 | // two hit times. So, all timestamps contain the overall IceEvent | |
74 | // timestamp as a basis. This means that time differences can be | |
75 | // obtained via the AliTimestamp facilities (supporting upto picosecond | |
76 | // precision when available). | |
77 | // The TE direction is given by the relative position of the two hits. | |
78 | // | |
79 | // 2) Each TE will obtain so called associated hits. | |
80 | // A hit is associated to a TE when it fulfills both the conditions | |
81 | // | |
82 | // -30 < tres < 300 ns | |
83 | // dhit < 25*(tres+30)^(1/4) meter | |
84 | // | |
85 | // tres : time residual | |
86 | // Difference between the observed hit time and the time expected | |
87 | // for a direct photon hit. | |
88 | // dhit : Distance between the hit and the TE | |
89 | // | |
90 | // 3) Construction of track candidates (TC's). | |
91 | // These are TE's that fulfill both the conditions (see eq. (21) in the NIM article) | |
92 | // | |
93 | // sigmal >= 20 meter | |
94 | // qtc >= 0.7*qtcmax | |
95 | // | |
96 | // qtc=min(nah,0.3*sigmal+7) | |
97 | // qtcmax=max(qtc) of all TE's with sigmal>=20 meter | |
98 | // | |
99 | // where we have used the observables : | |
100 | // | |
101 | // nah : Number of associated hits. | |
102 | // sigmal : rms variance of the distances between r0 and the projection | |
103 | // points on the track of the various associated hit positions. | |
104 | // | |
105 | // Note : The following additional quality selection as indicated | |
106 | // in the NIM article is not used anymore. | |
107 | // | |
108 | // nah >= 10 | |
109 | // | |
110 | // 4) The remaining track candidates are clustered into jets when their directions | |
111 | // are within a certain maximum opening angle. | |
112 | // In addition the distance between their r0's must be below a certain maximum | |
113 | // or the relative r0 direction must fall within a certain maximum opening angle | |
114 | // w.r.t. the jet-starting track candidate. | |
115 | // The latter criterion prevents clustering of (nearly) parallel track candidates | |
116 | // crossing the detector a very different locations (e.g. muon bundles). | |
117 | // The default maximum track opening angle is 15 degrees, but can be modified | |
118 | // via the SetTangmax memberfunction. | |
119 | // The remaining parameters related to the r0 criteria can be modified via | |
120 | // the SetRtdmax and SetRtangmax memberfunctions. | |
121 | // See the corresponding docs for defaults etc... | |
122 | // | |
123 | // The average of all the r0 and t0 values of the constituent TC's | |
124 | // of the jet will provide the r0 and t0 (i.e. reference point) of the jet. | |
125 | // | |
126 | // 5) The jets are merged when their directions are within a certain maximum | |
127 | // opening angle. Before the merging, all jets are ordered w.r.t. decreasing | |
128 | // number of associated hits. This guarantees that the jet with the largest | |
129 | // number of associated hits is the starting jet in the merging process. | |
130 | // In addition the distance between their r0's must be below a certain maximum | |
131 | // or the relative r0 direction must fall within a certain maximum opening angle | |
132 | // w.r.t. the starting jet. | |
133 | // The latter criterion prevents merging of (nearly) parallel tracks/jets | |
134 | // crossing the detector a very different locations (e.g. muon bundles). | |
135 | // The default maximum opening angle is half the TC maximum opening angle, | |
136 | // but can be modified via the SetJangmax memberfunction. | |
137 | // The remaining parameters related to the r0 criteria can be modified via | |
138 | // the SetRjdmax and SetRjangmax memberfunctions. | |
139 | // See the corresponding docs for defaults etc... | |
140 | // | |
141 | // Note : Setting the maximum jet opening angle to <=0 will prevent | |
142 | // the merging of jets. | |
143 | // | |
144 | // The average of all the r0 and t0 values of the merged jets will provide | |
145 | // the r0 and t0 (i.e. reference point) of the final jet. | |
146 | // | |
147 | // 6) The remaining jets are ordered w.r.t. decreasing number of associated hits. | |
148 | // Each remaining jet will provide the parameters (e.g. direction) | |
149 | // for a reconstructed track. | |
150 | // The track 3-momentum is set to the total jet 3-momentum, normalised | |
151 | // to 1 GeV. The mass and charge of the track are left 0. | |
152 | // The reference point data of the jet will provide the r0 and t0 | |
153 | // (i.e. reference point) of the track. | |
154 | // | |
155 | // All these tracks will be stored in the IceEvent structure with as | |
156 | // default "IceDwalkx" as the name of the track. | |
157 | // This track name identifier can be modified by the user via the | |
158 | // SetTrackName() memberfunction. This will allow unique identification | |
159 | // of tracks which are produced when re-processing existing data with | |
160 | // different criteria. | |
161 | // By default the charge of the produced tracks is set to 0, since | |
162 | // no distinction can be made between positive or negative tracks. | |
163 | // However, the user can define the track charge by invokation | |
164 | // of the memberfunction SetCharge(). | |
165 | // This facility may be used to distinguish tracks produced by the | |
166 | // various reconstruction algorithms in a (3D) colour display | |
167 | // (see the class AliHelix for further details). | |
168 | // | |
169 | // Note : In case the maximum jet opening angle was specified <0, | |
170 | // only the jet with the maximum number of associated hits will | |
171 | // appear as a reconstructed track in the IceEvent structure. | |
172 | // This will allow comparison with the standard Sieglinde | |
173 | // single track direct walk reconstruction results. | |
174 | // | |
175 | // For further details the user is referred to NIM A524 (2004) 169. | |
176 | // | |
177 | // Note : This algorithm works best on data which has been calibrated | |
178 | // (IceCalibrate), cross talk corrected (IceXtalk) and cleaned | |
179 | // from noise hits etc. (IceCleanHits). | |
180 | // | |
181 | //--- Author: Nick van Eijndhoven 07-oct-2005 Utrecht University | |
182 | //- Modified: NvE $Date$ Utrecht University | |
183 | /////////////////////////////////////////////////////////////////////////// | |
184 | ||
185 | #include "IceDwalkx.h" | |
186 | #include "Riostream.h" | |
187 | ||
188 | ClassImp(IceDwalkx) // Class implementation to enable ROOT I/O | |
189 | ||
190 | IceDwalkx::IceDwalkx(const char* name,const char* title) : TTask(name,title) | |
191 | { | |
192 | // Default constructor. | |
193 | // The various reconstruction parameters are initialised to the values | |
194 | // as mentioned in NIM A524 (2004) 179-180. | |
195 | // The newly introduced angular separation parameter for jet merging | |
196 | // is initialised as half the value of the angular separation parameter | |
197 | // for track candidate clustering. | |
198 | fDmin=50; | |
199 | fDtmarg=30; | |
200 | fTangmax=15; | |
201 | fRtangmax=fTangmax; | |
202 | fRtdmax=0; | |
203 | fJangmax=fTangmax/2.; | |
204 | fRjangmax=fRtangmax; | |
205 | fRjdmax=fDmin; | |
206 | fMaxmodA=999; | |
207 | fMinmodA=0; | |
208 | fMaxhitsA=1; | |
209 | fVgroup=1; | |
210 | fTrackname="IceDwalkx"; | |
211 | fCharge=0; | |
212 | } | |
213 | /////////////////////////////////////////////////////////////////////////// | |
214 | IceDwalkx::~IceDwalkx() | |
215 | { | |
216 | // Default destructor. | |
217 | } | |
218 | /////////////////////////////////////////////////////////////////////////// | |
219 | void IceDwalkx::SetDmin(Float_t d) | |
220 | { | |
221 | // Set minimum hit distance (in m) to form a track element. | |
222 | // In the constructor the default has been set to 50 meter, in accordance | |
223 | // to eq.(20) of NIM A524 (2004) 179. | |
224 | fDmin=d; | |
225 | } | |
226 | /////////////////////////////////////////////////////////////////////////// | |
227 | void IceDwalkx::SetDtmarg(Int_t dt) | |
228 | { | |
229 | // Set maximum hit time difference margin (in ns) for track elements. | |
230 | // In the constructor the default has been set to 30 ns, in accordance | |
231 | // to eq.(20) of NIM A524 (2004) 179. | |
232 | fDtmarg=dt; | |
233 | } | |
234 | /////////////////////////////////////////////////////////////////////////// | |
235 | void IceDwalkx::SetTangmax(Float_t ang) | |
236 | { | |
237 | // Set maximum angular separation (in deg) for track candidate clustering | |
238 | // into jets. | |
239 | // In the constructor the default has been set to 15 deg, in accordance | |
240 | // to NIM A524 (2004) 180. | |
241 | // | |
242 | // Note : This function also sets automatically the value of the maximum | |
243 | // angular separation for jet merging into 1 single track to ang/2. | |
244 | // In order to specify a different max. jet merging separation angle, | |
245 | // one has to invoke the memberfunction SetJangmax afterwards. | |
246 | ||
247 | fTangmax=ang; | |
248 | fJangmax=ang/2.; | |
249 | } | |
250 | /////////////////////////////////////////////////////////////////////////// | |
251 | void IceDwalkx::SetRtangmax(Float_t ang) | |
252 | { | |
253 | // Set maximum angular separation (in deg) for the relative direction of the | |
254 | // two r0's of two track candidates (w.r.t. the direction of the starting | |
255 | // track candidate) in the track clustering process. | |
256 | // In the constructor the default has been set to 15 deg, corresponding | |
257 | // to the default max. angular separation for track candidate clustering. | |
258 | // | |
259 | // Note : This function also sets automatically the value of the maximum | |
260 | // angular separation for the relative direction of the two r0's | |
261 | // of two jets (w.r.t. the direction of the starting jet) | |
262 | // in the jet merging process. | |
263 | // In order to specify a different value related to jet merging, | |
264 | // one has to invoke the memberfunction SetRjangmax afterwards. | |
265 | ||
266 | fRtangmax=ang; | |
267 | fRjangmax=ang; | |
268 | } | |
269 | /////////////////////////////////////////////////////////////////////////// | |
270 | void IceDwalkx::SetRtdmax(Float_t d) | |
271 | { | |
272 | // Set maximum distance (in m) of the two r0's of two track candidates | |
273 | // in the track clustering process. | |
274 | // This will allow clustering of tracks with very close r0's, of which | |
275 | // their relative direction may point in any direction. | |
276 | // In the constructor the default has been set 0 until further tuning | |
277 | // of this parameter has been achieved. | |
278 | // | |
279 | // Note : In case the distance between the two r0's exceeds the maximum, | |
280 | // the track candidates will still be clustered if the relative | |
281 | // direction of the two r0's falls within the maximum separation | |
282 | // angle w.r.t. the starting track direction. | |
283 | ||
284 | fRtdmax=d; | |
285 | } | |
286 | /////////////////////////////////////////////////////////////////////////// | |
287 | void IceDwalkx::SetJangmax(Float_t ang) | |
288 | { | |
289 | // Set angular separation (in deg) within which jets are merged into 1 | |
290 | // single track. | |
291 | // In the constructor the default has been set 7.5 deg, being half of the | |
292 | // value of the default track candidate clustering separation angle. | |
293 | // | |
294 | // Notes : | |
295 | // ------- | |
296 | // 1) Setting ang=0 will prevent jet merging. | |
297 | // Consequently, every jet will appear as a separate track in the | |
298 | // reconstruction result. | |
299 | // 2) Setting ang<0 will prevent jet merging. | |
300 | // In addition, only the jet with the maximum number of tracks will | |
301 | // appear as a track in the reconstruction result. | |
302 | // This situation resembles the standard Sieglinde direct walk processing | |
303 | // and as such can be used to perform comparison studies. | |
304 | ||
305 | fJangmax=ang; | |
306 | } | |
307 | /////////////////////////////////////////////////////////////////////////// | |
308 | void IceDwalkx::SetRjangmax(Float_t ang) | |
309 | { | |
310 | // Set maximum angular separation (in deg) for the relative direction of the | |
311 | // two r0's of two jets (w.r.t. the direction of the starting jet) | |
312 | // in the jet merging process. | |
313 | // In the constructor the default has been set to the corresponding value | |
314 | // of the same parameter related to the track clustering. | |
315 | ||
316 | fRjangmax=ang; | |
317 | } | |
318 | /////////////////////////////////////////////////////////////////////////// | |
319 | void IceDwalkx::SetRjdmax(Float_t d) | |
320 | { | |
321 | // Set maximum distance (in m) of the two r0's of two jets in the | |
322 | // jet merging process. | |
323 | // This will allow merging of jets with rather close r0's, of which | |
324 | // their relative direction may point in any direction. | |
325 | // In the constructor the default has been set to 50 meter, corresponding | |
326 | // to the value of the minimum hit distance to form a track element. | |
327 | // | |
328 | // Note : In case the distance between the two r0's exceeds the maximum, | |
329 | // the jets will still be merged if the relative direction of the | |
330 | // two r0's falls within the maximum separation angle w.r.t. the | |
331 | // starting jet direction. | |
332 | ||
333 | fRjdmax=d; | |
334 | } | |
335 | /////////////////////////////////////////////////////////////////////////// | |
336 | void IceDwalkx::SetMaxModA(Int_t nmax) | |
337 | { | |
338 | // Set the maximum number of good Amanda modules that may have fired | |
339 | // in order to process this event. | |
340 | // This allows suppression of processing (high-energy) cascade events | |
341 | // with this direct walk tracking to prevent waisting cpu time for cases | |
342 | // in which tracking doesn't make sense anyhow. | |
343 | // Furthermore it allows selection of low multiplicity events for processing. | |
344 | // By default the maximum number of Amanda modules is set to 999 in the ctor, | |
345 | // which implies no selection on maximum module multiplicity. | |
346 | // See also the memberfunction SetMinModA(). | |
347 | fMaxmodA=nmax; | |
348 | } | |
349 | /////////////////////////////////////////////////////////////////////////// | |
350 | void IceDwalkx::SetMinModA(Int_t nmin) | |
351 | { | |
352 | // Set the minimum number of good Amanda modules that must have fired | |
353 | // in order to process this event. | |
354 | // This allows selection of a minimal multiplicity for events to be processed. | |
355 | // By default the minimum number of Amanda modules is set to 0 in the ctor, | |
356 | // which implies no selection on minimum module multiplicity. | |
357 | // See also the memberfunction SetMaxModA(). | |
358 | fMinmodA=nmin; | |
359 | } | |
360 | /////////////////////////////////////////////////////////////////////////// | |
361 | void IceDwalkx::SetMaxHitsA(Int_t nmax) | |
362 | { | |
363 | // Set the maximum number of good hits per Amanda module to be processed. | |
364 | // | |
365 | // Special values : | |
366 | // nmax = 0 : No maximum limit set; all available good hits will be processed | |
367 | // < 0 : No hits will be processed | |
368 | // | |
369 | // In case the user selects a maximum number of good hits per module, all the | |
370 | // hits of each module will be ordered w.r.t. increasing hit time (LE). | |
371 | // This allows selection of processing e.g. only the first good hits etc... | |
372 | // By default the maximum number of good hits per Amanda modules is set to 1 | |
373 | // in the ctor, which implies just processing only the first good hit of | |
374 | // each Amanda OM. | |
375 | fMaxhitsA=nmax; | |
376 | } | |
377 | /////////////////////////////////////////////////////////////////////////// | |
378 | void IceDwalkx::SetVgroupUsage(Int_t flag) | |
379 | { | |
380 | // (De)activate the distinction between v_phase and v_group of the Cherenkov light. | |
381 | // | |
382 | // flag = 0 : No distinction between v_phase and v_group | |
383 | // = 1 : Separate treatment of v_phase and v_group | |
384 | // | |
385 | // By default the distinction between v_phase and v_group is activated | |
386 | // in the constructor of this class. | |
387 | fVgroup=flag; | |
388 | } | |
389 | /////////////////////////////////////////////////////////////////////////// | |
390 | void IceDwalkx::SetTrackName(TString s) | |
391 | { | |
392 | // Set (alternative) name identifier for the produced first guess tracks. | |
393 | // This allows unique identification of (newly) produced direct walk tracks | |
394 | // in case of re-processing of existing data with different criteria. | |
395 | // By default the produced first guess tracks have the name "IceDwalkx" | |
396 | // which is set in the constructor of this class. | |
397 | fTrackname=s; | |
398 | } | |
399 | /////////////////////////////////////////////////////////////////////////// | |
400 | void IceDwalkx::SetCharge(Float_t charge) | |
401 | { | |
402 | // Set user defined charge for the produced first guess tracks. | |
403 | // This allows identification of these tracks on color displays. | |
404 | // By default the produced first guess tracks have charge=0 | |
405 | // which is set in the constructor of this class. | |
406 | fCharge=charge; | |
407 | } | |
408 | /////////////////////////////////////////////////////////////////////////// | |
409 | void IceDwalkx::Exec(Option_t* opt) | |
410 | { | |
411 | // Implementation of the direct walk track reconstruction. | |
412 | ||
413 | TString name=opt; | |
414 | AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data())); | |
415 | ||
416 | if (!parent) return; | |
417 | ||
418 | IceEvent* evt=(IceEvent*)parent->GetObject("IceEvent"); | |
419 | if (!evt) return; | |
420 | ||
421 | // Enter the reco parameters as a device in the event | |
422 | AliSignal params; | |
423 | params.SetNameTitle("IceDwalkx","IceDwalkx reco parameters"); | |
424 | params.SetSlotName("Dmin",1); | |
425 | params.SetSlotName("Dtmarg",2); | |
426 | params.SetSlotName("Tangmax",3); | |
427 | params.SetSlotName("Rtangmax",4); | |
428 | params.SetSlotName("Rtdmax",5); | |
429 | params.SetSlotName("Jangmax",6); | |
430 | params.SetSlotName("Rjangmax",7); | |
431 | params.SetSlotName("Rjdmax",8); | |
432 | params.SetSlotName("MaxmodA",9); | |
433 | params.SetSlotName("MinmodA",10); | |
434 | params.SetSlotName("MaxhitsA",11); | |
435 | params.SetSlotName("Vgroup",12); | |
436 | ||
437 | params.SetSignal(fDmin,1); | |
438 | params.SetSignal(fDtmarg,2); | |
439 | params.SetSignal(fTangmax,3); | |
440 | params.SetSignal(fRtangmax,4); | |
441 | params.SetSignal(fRtdmax,5); | |
442 | params.SetSignal(fJangmax,6); | |
443 | params.SetSignal(fRjangmax,7); | |
444 | params.SetSignal(fRjdmax,8); | |
445 | params.SetSignal(fMaxmodA,9); | |
446 | params.SetSignal(fMinmodA,10); | |
447 | params.SetSignal(fMaxhitsA,11); | |
448 | params.SetSignal(fVgroup,12); | |
449 | ||
450 | evt->AddDevice(params); | |
451 | ||
452 | if (fMaxhitsA<0) return; | |
453 | ||
454 | // Fetch all fired Amanda OMs for this event | |
455 | TObjArray* aoms=evt->GetDevices("IceAOM"); | |
8a394508 | 456 | if (!aoms) return; |
25eefd00 | 457 | Int_t naoms=aoms->GetEntries(); |
458 | if (!naoms) return; | |
459 | ||
460 | // Check for the minimum and/or maximum number of good fired Amanda OMs | |
461 | Int_t ngood=0; | |
462 | for (Int_t iom=0; iom<naoms; iom++) | |
463 | { | |
464 | IceGOM* omx=(IceGOM*)aoms->At(iom); | |
465 | if (!omx) continue; | |
466 | if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue; | |
467 | ngood++; | |
468 | } | |
469 | if (ngood<fMinmodA || ngood>fMaxmodA) return; | |
470 | ||
471 | const Float_t pi=acos(-1.); | |
472 | const Float_t c=0.299792458; // Light speed in vacuum in meters per ns | |
473 | const Float_t npice=1.31768387; // Phase refractive index (c/v_phase) of ice | |
474 | const Float_t ngice=1.35075806; // Group refractive index (c/v_group) of ice | |
475 | const Float_t thetac=acos(1./npice); // Cherenkov angle (in radians) | |
476 | ||
477 | // Angular reduction of complement of thetac due to v_phase and v_group difference | |
478 | Float_t alphac=0; | |
479 | if (fVgroup) alphac=atan((1.-npice/ngice)/sqrt(npice*npice-1.)); | |
480 | ||
481 | // Storage of track elements. | |
482 | TObjArray tes; | |
483 | tes.SetOwner(); | |
484 | ||
485 | AliPosition r1; | |
486 | AliPosition r2; | |
487 | Ali3Vector r12; | |
488 | Ali3Vector rsum; | |
489 | AliPosition r0; | |
490 | TObjArray hits1; | |
491 | TObjArray hits2; | |
492 | Int_t nh1,nh2; | |
493 | AliSignal* sx1=0; | |
494 | AliSignal* sx2=0; | |
495 | Float_t dist=0; | |
496 | Float_t t1,t2,dt,t0; | |
497 | Float_t dtmax; | |
498 | TObjArray hits; | |
499 | TObjArray* ordered; | |
500 | ||
501 | // Check the hits of Amanda OM pairs for posible track elements. | |
502 | // Also all the good hits are stored in the meantime (to save CPU time) | |
503 | // for hit association with the various track elements lateron. | |
504 | AliTrack* te=0; | |
505 | for (Int_t i1=0; i1<naoms; i1++) // First OM of the pair | |
506 | { | |
507 | IceGOM* omx1=(IceGOM*)aoms->At(i1); | |
508 | if (!omx1) continue; | |
509 | if (omx1->GetDeadValue("LE")) continue; | |
510 | r1=omx1->GetPosition(); | |
511 | // Select all the good hits of this first OM | |
512 | hits1.Clear(); | |
513 | // Determine the max. number of hits to be processed for this OM | |
514 | ordered=0; | |
515 | if (fMaxhitsA>0 && omx1->GetNhits()>fMaxhitsA) ordered=omx1->SortHits("LE",1,0,7); | |
516 | nh1=0; | |
517 | for (Int_t j1=1; j1<=omx1->GetNhits(); j1++) | |
518 | { | |
519 | if (ordered) | |
520 | { | |
521 | if (nh1>=fMaxhitsA) break; | |
522 | sx1=(AliSignal*)ordered->At(j1-1); | |
523 | } | |
524 | else | |
525 | { | |
526 | sx1=omx1->GetHit(j1); | |
527 | } | |
528 | if (!sx1) continue; | |
529 | if (sx1->GetDeadValue("ADC") || sx1->GetDeadValue("LE") || sx1->GetDeadValue("TOT")) continue; | |
530 | hits1.Add(sx1); | |
531 | // Also store all good hits in the total hit array | |
532 | hits.Add(sx1); | |
533 | nh1++; | |
534 | } | |
535 | ||
536 | // No further pair to be formed with the last OM in the list | |
537 | if (i1==(naoms-1)) break; | |
538 | ||
539 | nh1=hits1.GetEntries(); | |
540 | if (!nh1) continue; | |
541 | ||
542 | for (Int_t i2=i1+1; i2<naoms; i2++) // Second OM of the pair | |
543 | { | |
544 | IceGOM* omx2=(IceGOM*)aoms->At(i2); | |
545 | if (!omx2) continue; | |
546 | if (omx2->GetDeadValue("LE")) continue; | |
547 | r2=omx2->GetPosition(); | |
548 | r12=r2-r1; | |
549 | dist=r12.GetNorm(); | |
550 | ||
551 | if (dist<=fDmin) continue; | |
552 | ||
553 | // Select all the good hits of this second OM | |
554 | hits2.Clear(); | |
555 | // Determine the max. number of hits to be processed for this OM | |
556 | ordered=0; | |
557 | if (fMaxhitsA>0 && omx2->GetNhits()>fMaxhitsA) ordered=omx2->SortHits("LE",1,0,7); | |
558 | nh2=0; | |
559 | for (Int_t j2=1; j2<=omx2->GetNhits(); j2++) | |
560 | { | |
561 | if (ordered) | |
562 | { | |
563 | if (nh2>=fMaxhitsA) break; | |
564 | sx2=(AliSignal*)ordered->At(j2-1); | |
565 | } | |
566 | else | |
567 | { | |
568 | sx2=omx2->GetHit(j2); | |
569 | } | |
570 | if (!sx2) continue; | |
571 | if (sx2->GetDeadValue("ADC") || sx2->GetDeadValue("LE") || sx2->GetDeadValue("TOT")) continue; | |
572 | hits2.Add(sx2); | |
573 | nh2++; | |
574 | } | |
575 | ||
576 | nh2=hits2.GetEntries(); | |
577 | if (!nh2) continue; | |
578 | ||
579 | // Position r0 in between the two OMs and normalised relative direction r12 | |
580 | rsum=(r1+r2)/2.; | |
581 | r0.SetPosition((Ali3Vector&)rsum); | |
582 | r12/=dist; | |
583 | ||
584 | // Check all hit pair combinations of these two OMs for possible track elements | |
585 | dtmax=dist/c+float(fDtmarg); | |
586 | for (Int_t ih1=0; ih1<nh1; ih1++) // Hits of first OM | |
587 | { | |
588 | sx1=(AliSignal*)hits1.At(ih1); | |
589 | if (!sx1) continue; | |
590 | for (Int_t ih2=0; ih2<nh2; ih2++) // Hits of second OM | |
591 | { | |
592 | sx2=(AliSignal*)hits2.At(ih2); | |
593 | if (!sx2) continue; | |
594 | t1=sx1->GetSignal("LE",7); | |
595 | t2=sx2->GetSignal("LE",7); | |
596 | dt=t2-t1; | |
597 | t0=(t1+t2)/2.; | |
598 | ||
599 | if (fabs(dt)>=dtmax) continue; | |
600 | ||
601 | te=new AliTrack(); | |
602 | tes.Add(te); | |
603 | if (dt<0) r12*=-1.; | |
604 | r0.SetTimestamp((AliTimestamp&)*evt); | |
605 | AliTimestamp* tsx=r0.GetTimestamp(); | |
606 | tsx->Add(0,0,(int)t0); | |
607 | te->SetReferencePoint(r0); | |
608 | te->Set3Momentum(r12); | |
609 | } | |
610 | } | |
611 | } // end of loop over the second OM of the pair | |
612 | } // end of loop over first OM of the pair | |
613 | ||
614 | // Association of hits to the various track elements | |
615 | // For the time being all track elements will be treated, | |
616 | // but in a later stage one could select only the TE's of a certain | |
617 | // 3 ns margin slot in the TE map to save CPU time. | |
618 | Int_t nte=tes.GetEntries(); | |
619 | Int_t nh=hits.GetEntries(); | |
620 | Float_t d=0; | |
621 | Ali3Vector p; | |
622 | Float_t tgeo,tres; | |
623 | AliSample levers; // Statistics of the assoc. hit lever arms | |
624 | AliSignal fit; // Storage of Q value etc... for each track candidate | |
625 | fit.SetSlotName("QTC",1); | |
626 | fit.SetSlotName("SIGMAL",2); | |
627 | Float_t qtc=0,qmax=0; | |
628 | Int_t nah; // Number of associated hits for a certain TE | |
629 | Float_t sigmal; // The mean lever arm of the various associated hits | |
630 | Float_t hproj; // Projected hit position on the track w.r.t. r0 | |
631 | for (Int_t jte=0; jte<nte; jte++) | |
632 | { | |
633 | te=(AliTrack*)tes.At(jte); | |
634 | if (!te) continue; | |
635 | AliPosition* tr0=te->GetReferencePoint(); | |
636 | AliTimestamp* tt0=tr0->GetTimestamp(); | |
637 | t0=evt->GetDifference(tt0,"ns"); | |
638 | p=te->Get3Momentum(); | |
639 | levers.Reset(); | |
640 | for (Int_t jh=0; jh<nh; jh++) | |
641 | { | |
642 | sx1=(AliSignal*)hits.At(jh); | |
643 | if (!sx1) continue; | |
644 | IceGOM* omx=(IceGOM*)sx1->GetDevice(); | |
645 | if (!omx) continue; | |
646 | r1=omx->GetPosition(); | |
647 | d=te->GetDistance(r1); | |
648 | r12=r1-(*tr0); | |
649 | hproj=p.Dot(r12); | |
650 | dist=hproj+d/tan(pi/2.-thetac-alphac); | |
651 | tgeo=t0+dist/c; | |
652 | t1=sx1->GetSignal("LE",7); | |
653 | tres=t1-tgeo; | |
654 | ||
655 | if (tres<-30 || tres>300 || d>25.*pow(tres+30.,0.25)) continue; | |
656 | ||
657 | // Associate this hit to the TE | |
658 | te->AddSignal(*sx1); | |
659 | levers.Enter(fabs(hproj)); | |
660 | } | |
661 | ||
662 | // Determine the Q quality of the various TE's. | |
663 | // Good quality TE's will be called track candidates (TC's) | |
664 | nah=te->GetNsignals(); | |
665 | sigmal=levers.GetSigma(1); | |
666 | qtc=0.3*sigmal+7.; | |
667 | if (qtc>nah) qtc=nah; | |
668 | if (sigmal<20) qtc=-1; // Reject TE's with sigmal<20 meter | |
669 | fit.SetSignal(qtc,"QTC"); | |
670 | fit.SetSignal(sigmal,"SIGMAL"); | |
671 | te->SetFitDetails(fit); | |
672 | if (qtc>qmax) qmax=qtc; | |
673 | } | |
674 | ||
675 | // Perform selection on Q value in case of multiple track candidates | |
676 | for (Int_t jtc=0; jtc<nte; jtc++) | |
677 | { | |
678 | te=(AliTrack*)tes.At(jtc); | |
679 | if (!te) continue; | |
680 | nah=te->GetNsignals(); | |
681 | sx1=(AliSignal*)te->GetFitDetails(); | |
682 | qtc=-1; | |
683 | sigmal=0; | |
684 | if (sx1) | |
685 | { | |
686 | qtc=sx1->GetSignal("QTC"); | |
687 | sigmal=sx1->GetSignal("SIGMAL"); | |
688 | } | |
689 | if (qtc<0.7*qmax) | |
690 | { | |
691 | tes.RemoveAt(jtc); | |
692 | delete te; | |
693 | } | |
694 | } | |
695 | tes.Compress(); | |
696 | nte=tes.GetEntries(); | |
697 | ||
698 | if (!nte) return; | |
699 | ||
700 | // Cluster track candidates within a certain opening angle into jets. | |
701 | // Also the relative direction of the both r0's of the track candidates | |
702 | // should be within a certain opening angle w.r.t. the starting track direction, | |
703 | // unless the distance between the two r0's is below a certain maximum. | |
704 | // The latter prevents clustering of (nearly) parallel track candidates | |
705 | // crossing the detector a very different locations (e.g. muon bundles). | |
706 | // The average r0 and t0 of the constituent tracks will be taken as the | |
707 | // jet reference point. | |
708 | TObjArray jets; | |
709 | jets.SetOwner(); | |
710 | AliTrack* te2=0; | |
711 | Float_t ang=0; | |
712 | AliSample pos; | |
713 | AliSample time; | |
714 | Float_t vec[3],err[3]; | |
715 | for (Int_t jtc1=0; jtc1<nte; jtc1++) | |
716 | { | |
717 | te=(AliTrack*)tes.At(jtc1); | |
718 | if (!te) continue; | |
719 | AliPosition* x1=te->GetReferencePoint(); | |
720 | if (!x1) continue; | |
721 | AliTimestamp* ts1=x1->GetTimestamp(); | |
722 | if (!ts1) continue; | |
723 | AliJet* jx=new AliJet(); | |
724 | jx->AddTrack(te); | |
725 | pos.Reset(); | |
726 | time.Reset(); | |
727 | x1->GetPosition(vec,"car"); | |
728 | pos.Enter(vec[0],vec[1],vec[2]); | |
729 | t0=evt->GetDifference(ts1,"ns"); | |
730 | time.Enter(t0); | |
731 | for (Int_t jtc2=0; jtc2<nte; jtc2++) | |
732 | { | |
733 | if (jtc2==jtc1) continue; | |
734 | te2=(AliTrack*)tes.At(jtc2); | |
735 | if (!te2) continue; | |
736 | ang=te->GetOpeningAngle(*te2,"deg"); | |
737 | if (ang<=fTangmax) | |
738 | { | |
739 | AliPosition* x2=te2->GetReferencePoint(); | |
740 | if (!x2) continue; | |
741 | AliTimestamp* ts2=x2->GetTimestamp(); | |
742 | if (!ts2) continue; | |
743 | dist=x1->GetDistance(x2); | |
744 | dt=ts1->GetDifference(ts2,"ns"); | |
745 | if (dist>0) | |
746 | { | |
747 | r12=(*x2)-(*x1); | |
748 | if (dt<0) r12*=-1.; | |
749 | ang=te->GetOpeningAngle(r12,"deg"); | |
750 | if (ang<=fRtangmax || dist<fRtdmax) | |
751 | { | |
752 | x2->GetPosition(vec,"car"); | |
753 | pos.Enter(vec[0],vec[1],vec[2]); | |
754 | t0=evt->GetDifference(ts2,"ns"); | |
755 | time.Enter(t0); | |
756 | jx->AddTrack(te2); | |
757 | } | |
758 | } | |
759 | } | |
760 | } | |
761 | ||
762 | // Set the reference point data for this jet | |
763 | for (Int_t j=1; j<=3; j++) | |
764 | { | |
765 | vec[j-1]=pos.GetMean(j); | |
766 | err[j-1]=pos.GetSigma(j); | |
767 | } | |
768 | r0.SetPosition(vec,"car"); | |
769 | r0.SetPositionErrors(err,"car"); | |
770 | r0.SetTimestamp((AliTimestamp&)*evt); | |
771 | AliTimestamp* jt0=r0.GetTimestamp(); | |
772 | t0=time.GetMean(1); | |
773 | jt0->Add(0,0,(int)t0); | |
774 | jx->SetReferencePoint(r0); | |
775 | ||
776 | // Store this jet for further processing if ntracks>1 | |
777 | if (jx->GetNtracks() > 1 || fTangmax<=0 || fRtangmax<=0) | |
778 | { | |
779 | jets.Add(jx); | |
780 | } | |
781 | else // Only keep single-track jets which have qtc=qmax | |
782 | { | |
783 | sx1=(AliSignal*)te->GetFitDetails(); | |
784 | qtc=-1; | |
785 | if (sx1) qtc=sx1->GetSignal("QTC"); | |
786 | if (qtc>=(qmax-1.e-10)) | |
787 | { | |
788 | jets.Add(jx); | |
789 | } | |
790 | else | |
791 | { | |
792 | delete jx; | |
793 | } | |
794 | } | |
795 | } | |
796 | Int_t njets=jets.GetEntries(); | |
797 | ||
798 | if (!njets) return; | |
799 | ||
800 | // Order the jets w.r.t. decreasing number of associated hits | |
801 | ordered=evt->SortJets(-12,&jets); | |
802 | TObjArray jets2(*ordered); | |
803 | ||
804 | // Merge jets within a certain opening to provide the final track(s). | |
805 | AliJet* jx1=0; | |
806 | AliJet* jx2=0; | |
807 | Float_t edist=0; | |
808 | if (fJangmax>=0) | |
809 | { | |
810 | for (Int_t jet1=0; jet1<njets; jet1++) | |
811 | { | |
812 | jx1=(AliJet*)jets2.At(jet1); | |
813 | if (!jx1) continue; | |
814 | AliPosition* x1=jx1->GetReferencePoint(); | |
815 | if (!x1) continue; | |
816 | AliTimestamp* ts1=x1->GetTimestamp(); | |
817 | if (!ts1) continue; | |
818 | pos.Reset(); | |
819 | time.Reset(); | |
820 | x1->GetPosition(vec,"car"); | |
821 | pos.Enter(vec[0],vec[1],vec[2]); | |
822 | t0=evt->GetDifference(ts1,"ns"); | |
823 | time.Enter(t0); | |
824 | for (Int_t jet2=jet1+1; jet2<njets; jet2++) | |
825 | { | |
826 | jx2=(AliJet*)jets2.At(jet2); | |
827 | if (!jx2) continue; | |
828 | AliPosition* x2=jx2->GetReferencePoint(); | |
829 | if (!x2) continue; | |
830 | AliTimestamp* ts2=x2->GetTimestamp(); | |
831 | if (!ts2) continue; | |
832 | ang=jx1->GetOpeningAngle(*jx2,"deg"); | |
833 | if (ang<=fJangmax) | |
834 | { | |
835 | dist=x1->GetDistance(x2); | |
836 | edist=x1->GetResultError(); | |
837 | dt=ts1->GetDifference(ts2,"ns"); | |
838 | r12=(*x2)-(*x1); | |
839 | if (dt<0) r12*=-1.; | |
840 | ang=jx1->GetOpeningAngle(r12,"deg"); | |
841 | if (ang<=fRjangmax || dist<=(fRjdmax+edist)) | |
842 | { | |
843 | x2->GetPosition(vec,"car"); | |
844 | pos.Enter(vec[0],vec[1],vec[2]); | |
845 | t0=evt->GetDifference(ts2,"ns"); | |
846 | time.Enter(t0); | |
847 | for (Int_t jtk=1; jtk<=jx2->GetNtracks(); jtk++) | |
848 | { | |
849 | te=jx2->GetTrack(jtk); | |
850 | if (!te) continue; | |
851 | jx1->AddTrack(te); | |
852 | } | |
853 | jets2.RemoveAt(jet2); | |
854 | } | |
855 | } | |
856 | } | |
857 | // Set the reference point data for this jet | |
858 | for (Int_t k=1; k<=3; k++) | |
859 | { | |
860 | vec[k-1]=pos.GetMean(k); | |
861 | err[k-1]=pos.GetSigma(k); | |
862 | } | |
863 | r0.SetPosition(vec,"car"); | |
864 | r0.SetPositionErrors(err,"car"); | |
865 | r0.SetTimestamp((AliTimestamp&)*evt); | |
866 | AliTimestamp* jt0=r0.GetTimestamp(); | |
867 | t0=time.GetMean(1); | |
868 | jt0->Add(0,0,(int)t0); | |
869 | jx1->SetReferencePoint(r0); | |
870 | } | |
871 | jets2.Compress(); | |
872 | njets=jets2.GetEntries(); | |
873 | } | |
874 | ||
875 | // Order the merged jets w.r.t. decreasing number of associated hits | |
876 | ordered=evt->SortJets(-12,&jets2); | |
877 | TObjArray jets3(*ordered); | |
878 | ||
879 | // Store every jet as a reconstructed track in the event structure. | |
880 | // The jet 3-momentum (normalised to 1) and reference point | |
881 | // (i.e.the average r0 and t0 of the constituent tracks) will make up | |
882 | // the final track parameters. | |
883 | // All the associated hits of all the constituent tracks of the jet | |
884 | // will be associated to the final track. | |
885 | // In case the jet angular separation was set <0, only the jet with | |
886 | // the maximum number of tracks (i.e. the first one in the array) | |
887 | // will be used to form a track. This will allow comparison with | |
888 | // the standard Sieglinde processing. | |
889 | AliTrack t; | |
890 | t.SetNameTitle(fTrackname.Data(),"IceDwalkx direct walk track"); | |
891 | t.SetCharge(fCharge); | |
892 | for (Int_t jet=0; jet<njets; jet++) | |
893 | { | |
894 | AliJet* jx=(AliJet*)jets3.At(jet); | |
895 | if (!jx) continue; | |
896 | AliPosition* ref=jx->GetReferencePoint(); | |
897 | if (!ref) continue; | |
898 | evt->AddTrack(t); | |
899 | AliTrack* trk=evt->GetTrack(evt->GetNtracks()); | |
900 | if (!trk) continue; | |
901 | trk->SetId(evt->GetNtracks(1)+1); | |
902 | p=jx->Get3Momentum(); | |
903 | p/=p.GetNorm(); | |
904 | trk->Set3Momentum(p); | |
905 | trk->SetReferencePoint(*ref); | |
906 | AliTimestamp* tt0=ref->GetTimestamp(); | |
907 | if (tt0) trk->SetTimestamp(*tt0); | |
908 | for (Int_t jt=1; jt<=jx->GetNtracks(); jt++) | |
909 | { | |
910 | AliTrack* tx=jx->GetTrack(jt); | |
911 | if (!tx) continue; | |
912 | for (Int_t is=1; is<=tx->GetNsignals(); is++) | |
913 | { | |
914 | sx1=tx->GetSignal(is); | |
915 | if (sx1) sx1->AddTrack(*trk); | |
916 | } | |
917 | } | |
918 | ||
919 | // Only take the jet with the maximum number of associated hits | |
920 | // (i.e. the first jet in the list) when the user had selected | |
921 | // this reconstruction mode. | |
922 | if (fJangmax<0) break; | |
923 | } | |
924 | } | |
925 | /////////////////////////////////////////////////////////////////////////// |