]>
Commit | Line | Data |
---|---|---|
67338e65 | 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 IceDwalk | |
20 | // TTask derived class to perform direct walk track reconstruction. | |
caa58e1a | 21 | // |
22 | // In case an event has been rejected by an AliEventSelector (based) processor, | |
23 | // this task (and its sub-tasks) is not executed. | |
24 | // | |
67338e65 | 25 | // The procedure is based on the method described in the Amanda publication |
26 | // in Nuclear Instruments and Methods A524 (2004) 179-180. | |
27 | // However, the Amanda method has been extended with the intention to | |
28 | // take also multiple (muon) tracks within 1 event into account. | |
29 | // This will not only provide a means to reconstruct muon bundles and | |
30 | // multiple track events in IceCube, but will also allow to reduce the | |
31 | // background of faked upgoing muons as a result of multiple downgoing | |
d6860cf1 | 32 | // muons hitting the top and bottom parts of the detector. |
25eefd00 | 33 | // A further extension of the original Amanda method is the separate treatment |
34 | // of the phase and group velocities as introduced in collaboration with | |
35 | // George Japaridze (Clark Atlanta University, USA) which will provide more | |
36 | // accurate time residuals due to the different velocities of the Cerenkov | |
37 | // wave front (v_phase) and the actually detected photons (v_group). | |
38 | // This distinction between v_phase and v_group can be (de)activated via the | |
39 | // memberfunction SetVgroupUsage(). By default the distinction between v_phase | |
40 | // and v_group is activated in the constructor of this class. | |
d6860cf1 | 41 | // To prevent waisting CPU time in trying to reconstruct (high-energy) cascade |
42 | // events, or to select specifically reconstruction of low multiplicity events, | |
6f94f699 | 43 | // the user may invoke the memberfunctions SetMaxModA() and SetMinModA(). |
44 | // This allows selection of events for processing with a certain maximum and/or | |
45 | // minimum number of good Amanda OMs firing. | |
46 | // By default the minimum and maximum are set to 0 and 999, respectively, | |
25eefd00 | 47 | // in the constructor, which implies no multiplicity selection. |
48 | // The maximum number of good hits per Amanda OM to be used for the reconstruction | |
35e721bf | 49 | // can be specified via the memberfunction SetMaxHitsA(). |
25eefd00 | 50 | // By default only the first good hit of each Amanda OM is used. |
51 | // Note that when all the good hits of an OM are used, this may lead to large | |
35e721bf | 52 | // processing time in case many noise and/or afterpulse signals are not |
53 | // recognised by the hit cleaning procedure. | |
25eefd00 | 54 | // |
55 | // Information about the actual parameter settings can be found in the event | |
56 | // structure itself via the device named "IceDwalk". | |
57 | // | |
67338e65 | 58 | // The various reconstruction steps are summarised as follows : |
59 | // | |
60 | // 1) Construction of track elements (TE's). | |
61 | // A track element is a straight line connecting two hits that | |
62 | // appeared at some minimum distance d and within some maximum | |
25eefd00 | 63 | // time difference dt, according to eq. (20) of the NIM article. |
64 | // The default value for d is 75 meter, but this can be modified | |
65 | // via the memberfunction SetDmin(). | |
66 | // By default dt=(hit distance)/c but an additional time margin | |
67 | // may be specified via the memberfunction SetDtmarg(). | |
67338e65 | 68 | // The reference point r0 of the TE is taken as the center between |
69 | // the two hit positions and the TE timestamp t0 at the position r0 | |
70 | // is taken as the IceEvent timestamp increased by the average of the | |
71 | // two hit times. So, all timestamps contain the overall IceEvent | |
72 | // timestamp as a basis. This means that time differences can be | |
73 | // obtained via the AliTimestamp facilities (supporting upto picosecond | |
74 | // precision when available). | |
75 | // The TE direction is given by the relative position of the two hits. | |
76 | // | |
77 | // 2) Each TE will obtain so called associated hits. | |
78 | // A hit is associated to a TE when it fulfills both the conditions | |
79 | // | |
80 | // -30 < tres < 300 ns | |
25eefd00 | 81 | // dhit/lambda < F |
82 | // | |
83 | // tres : time residual | |
84 | // Difference between the observed hit time and the time expected | |
85 | // for a direct photon hit. | |
bfde1ea0 | 86 | // dhit : Distance traveled by the cherenkov photon from the track to the hit position |
25eefd00 | 87 | // lambda : Photon scattering length in ice |
67338e65 | 88 | // |
bfde1ea0 | 89 | // By default F is set to 3.07126, but this can be modified via the memberfunction |
25eefd00 | 90 | // SetMaxDhit(). |
67338e65 | 91 | // |
92 | // 3) Construction of track candidates (TC's). | |
25eefd00 | 93 | // These are TE's that fulfill both the conditions |
94 | // | |
95 | // nah >= 1 | |
96 | // qtc >= 0.8*qtcmax | |
97 | // | |
98 | // where we have defined : | |
413d0114 | 99 | // |
25eefd00 | 100 | // nah : Number of associated hits for the specific TE. |
101 | // qtc : The track quality number (see hereafter). | |
102 | // qtcmax : Maximum quality number encountered for the TE's. | |
413d0114 | 103 | // |
25eefd00 | 104 | // The track quality number qtc is defined as follows : |
413d0114 | 105 | // |
25eefd00 | 106 | // qtc=nah*(term1+term2)-term3-term4-term5 |
67338e65 | 107 | // |
25eefd00 | 108 | // here we have defined : |
67338e65 | 109 | // |
25eefd00 | 110 | // term1=2*spread/span |
111 | // term2=2*spreadL/spanL | |
112 | // term3=|spread-expspread|/spread | |
113 | // term4=|spreadL-expspreadL|/spreadL | |
114 | // term5=|medianT|/spreadT | |
67338e65 | 115 | // |
25eefd00 | 116 | // The central observables here are the projected positions X on the track |
117 | // of the various associated hits w.r.t. the track reference point r0. | |
118 | // Note that X can be negative as well as positive. | |
119 | // Therefore we also introduce XL=|X|. | |
120 | // | |
121 | // span : max(X)-min(X) | |
122 | // spanL : max(XL)-min(XL) | |
123 | // Xmedian : median of X | |
124 | // XmedianL : median of XL | |
125 | // spread : < |X-Xmedian| > | |
126 | // spreadL : < |XL-XmedianL| > | |
127 | // expspread : expected spread in X for a flat distribution of nah hits over span | |
128 | // expspreadL : expected spread in XL for a flat distribution of nah hits over spanL | |
129 | // medianT : median of tres | |
130 | // spreadT : < |tres-medianT| > | |
131 | // | |
132 | // However, if |Xmedian| > span/2 we set qtc=0 in order to always require | |
133 | // projected hits to appear on both sides of r0 on the track. | |
134 | // | |
135 | // Note : The qtc quality number is used to define the norm of the momentum | |
136 | // of the track candidate. As such it serves as a weight for the jet | |
137 | // momentum (direction) after clustering of the TC's and lateron | |
138 | // merging of the jets (see hereafter). | |
67338e65 | 139 | // |
413d0114 | 140 | // 4) The remaining track candidates are clustered into jets when their directions |
141 | // are within a certain maximum opening angle. | |
25eefd00 | 142 | // In addition a track candidate must within a certain maximum distance |
143 | // of the jet starting TC in order to get clustered. | |
413d0114 | 144 | // The latter criterion prevents clustering of (nearly) parallel track candidates |
145 | // crossing the detector a very different locations (e.g. muon bundles). | |
146 | // The default maximum track opening angle is 15 degrees, but can be modified | |
147 | // via the SetTangmax memberfunction. | |
25eefd00 | 148 | // The default maximum track distance is 20 meters, but can be modified |
149 | // via the SetTdistmax memberfunction. This memberfunction also allows to | |
150 | // specify whether the distance is determined within the detector volume or not. | |
67338e65 | 151 | // |
413d0114 | 152 | // The average of all the r0 and t0 values of the constituent TC's |
153 | // of the jet will provide the r0 and t0 (i.e. reference point) of the jet. | |
67338e65 | 154 | // |
25eefd00 | 155 | // The jet total momentum consists of the vector sum of the momenta of the |
156 | // constituent TC's. This implies that the qtc quality numbers of the various | |
157 | // TC's define a weight for each track in the construction of the jet direction. | |
158 | // In addition it means that the total jet momentum represents the sum of the | |
159 | // qtc quality numbers of the constituent TC's weighted by the opening angles | |
160 | // between the various TC's. | |
161 | // As such each jet is given an absolute quality number defined as : | |
162 | // | |
163 | // qtcjet=|jet momentum|/ntracks | |
164 | // | |
165 | // This jet quality number is refined on basis of the number of hits | |
166 | // associated to the jet as : | |
167 | // | |
168 | // qtcjet=qtcjet+0.2*(nah-nahmax) | |
169 | // | |
170 | // where we have defined : | |
171 | // | |
172 | // nah : Number of associated hits for the specific jet. | |
173 | // nahmax : Maximum number of associated hits encountered for the jets. | |
174 | // | |
175 | // This qtcjet value is then used to order the various jets w.r.t. | |
176 | // decreasing qtcjet quality number. | |
177 | // | |
178 | // Note : The qtcjet value is stored as "energy" of the jet, such that | |
179 | // it is always available for each jet and can also be used for | |
180 | // ordering the jets according to this value using the generic | |
181 | // AliEvent::SortJets() facility. | |
182 | // | |
183 | // 5) The jets (after having been ordered w.r.t. decreasing qtcjet value) | |
184 | // are merged when their directions are within a certain maximum opening angle. | |
185 | // In addition a jet must within a certain maximum distance of the starting jet | |
186 | // in order to get merged. | |
413d0114 | 187 | // The latter criterion prevents merging of (nearly) parallel tracks/jets |
188 | // crossing the detector a very different locations (e.g. muon bundles). | |
25eefd00 | 189 | // The jet ordering before the merging process is essential, since the starting jet |
190 | // will "eat up" the jets that will be merged into it. | |
191 | // The jet ordering ensures that the jet with the highest quality number will | |
192 | // always initiate the merging process. | |
67338e65 | 193 | // The default maximum opening angle is half the TC maximum opening angle, |
25eefd00 | 194 | // but can be modified via the SetJangmax memberfunction. This memberfunction |
195 | // also allows to specify whether jet merging will be performed iteratively or not. | |
196 | // In case iteration has been activated, the jet ordering is performed after each | |
197 | // iteration step. This has to be done because since the quality numbers of the | |
198 | // resulting merged jets have been automatically updated in the merging process. | |
199 | // | |
200 | // The default maximum jet distance is 30 meters, but can be modified | |
201 | // via the SetJdistmax memberfunction. This memberfunction also allows to | |
202 | // specify whether the distance is determined within the detector volume or not. | |
413d0114 | 203 | // |
67338e65 | 204 | // Note : Setting the maximum jet opening angle to <=0 will prevent |
205 | // the merging of jets. | |
206 | // | |
413d0114 | 207 | // The average of all the r0 and t0 values of the merged jets will provide |
208 | // the r0 and t0 (i.e. reference point) of the final jet. | |
209 | // | |
25eefd00 | 210 | // 6) The remaining (merged) jets are ordered w.r.t. decreasing jet quality number. |
211 | // As such the jet with the highest quality number will be the first one | |
212 | // in the list, which will result in the fact that the final tracks are also | |
213 | // ordered w.r.t. decreasing quality number, as outlined hereafter. | |
67338e65 | 214 | // Each remaining jet will provide the parameters (e.g. direction) |
215 | // for a reconstructed track. | |
216 | // The track 3-momentum is set to the total jet 3-momentum, normalised | |
217 | // to 1 GeV. The mass and charge of the track are left 0. | |
413d0114 | 218 | // The reference point data of the jet will provide the r0 and t0 |
219 | // (i.e. reference point) of the track. | |
d6860cf1 | 220 | // |
221 | // All these tracks will be stored in the IceEvent structure with as | |
222 | // default "IceDwalk" as the name of the track. | |
223 | // This track name identifier can be modified by the user via the | |
224 | // SetTrackName() memberfunction. This will allow unique identification | |
225 | // of tracks which are produced when re-processing existing data with | |
226 | // different criteria. | |
30672a96 | 227 | // By default the charge of the produced tracks is set to 0, since |
228 | // no distinction can be made between positive or negative tracks. | |
229 | // However, the user can define the track charge by invokation | |
230 | // of the memberfunction SetCharge(). | |
231 | // This facility may be used to distinguish tracks produced by the | |
232 | // various reconstruction algorithms in a (3D) colour display | |
233 | // (see the class AliHelix for further details). | |
413d0114 | 234 | // |
67338e65 | 235 | // Note : In case the maximum jet opening angle was specified <0, |
25eefd00 | 236 | // only the jet with the highest quality number will appear |
67338e65 | 237 | // as a reconstructed track in the IceEvent structure. |
238 | // This will allow comparison with the standard Sieglinde | |
239 | // single track direct walk reconstruction results. | |
240 | // | |
241 | // For further details the user is referred to NIM A524 (2004) 169. | |
242 | // | |
243 | // Note : This algorithm works best on data which has been calibrated | |
244 | // (IceCalibrate), cross talk corrected (IceXtalk) and cleaned | |
245 | // from noise hits etc. (IceCleanHits). | |
246 | // | |
247 | //--- Author: Nick van Eijndhoven 07-oct-2005 Utrecht University | |
248 | //- Modified: NvE $Date$ Utrecht University | |
249 | /////////////////////////////////////////////////////////////////////////// | |
250 | ||
251 | #include "IceDwalk.h" | |
252 | #include "Riostream.h" | |
253 | ||
254 | ClassImp(IceDwalk) // Class implementation to enable ROOT I/O | |
255 | ||
256 | IceDwalk::IceDwalk(const char* name,const char* title) : TTask(name,title) | |
257 | { | |
258 | // Default constructor. | |
259 | // The various reconstruction parameters are initialised to the values | |
260 | // as mentioned in NIM A524 (2004) 179-180. | |
261 | // The newly introduced angular separation parameter for jet merging | |
262 | // is initialised as half the value of the angular separation parameter | |
263 | // for track candidate clustering. | |
bfde1ea0 | 264 | fEvt=0; |
25eefd00 | 265 | fDmin=75; |
266 | fDtmarg=0; | |
bfde1ea0 | 267 | fMaxdhit=3.07126; |
413d0114 | 268 | fTangmax=15; |
25eefd00 | 269 | fTdistmax=20; |
270 | fTinvol=1; | |
413d0114 | 271 | fJangmax=fTangmax/2.; |
25eefd00 | 272 | fJiterate=1; |
273 | fJdistmax=30; | |
274 | fJinvol=1; | |
d6860cf1 | 275 | fMaxmodA=999; |
6f94f699 | 276 | fMinmodA=0; |
25eefd00 | 277 | fMaxhitsA=1; |
278 | fVgroup=1; | |
bfde1ea0 | 279 | fTrackname=""; |
30672a96 | 280 | fCharge=0; |
67338e65 | 281 | } |
282 | /////////////////////////////////////////////////////////////////////////// | |
283 | IceDwalk::~IceDwalk() | |
284 | { | |
285 | // Default destructor. | |
286 | } | |
287 | /////////////////////////////////////////////////////////////////////////// | |
288 | void IceDwalk::SetDmin(Float_t d) | |
289 | { | |
290 | // Set minimum hit distance (in m) to form a track element. | |
25eefd00 | 291 | // In the constructor the default has been set to 75 meter. |
67338e65 | 292 | fDmin=d; |
293 | } | |
294 | /////////////////////////////////////////////////////////////////////////// | |
295 | void IceDwalk::SetDtmarg(Int_t dt) | |
296 | { | |
297 | // Set maximum hit time difference margin (in ns) for track elements. | |
25eefd00 | 298 | // In the constructor the default has been set to 0 ns. |
67338e65 | 299 | fDtmarg=dt; |
300 | } | |
301 | /////////////////////////////////////////////////////////////////////////// | |
25eefd00 | 302 | void IceDwalk::SetMaxDhit(Float_t d) |
303 | { | |
304 | // Set maximum distance (in scattering length) for a hit to get associated. | |
305 | // In the constructor the default has been set to 2 lambda_scat. | |
306 | fMaxdhit=d; | |
307 | } | |
308 | /////////////////////////////////////////////////////////////////////////// | |
413d0114 | 309 | void IceDwalk::SetTangmax(Float_t ang) |
67338e65 | 310 | { |
413d0114 | 311 | // Set maximum angular separation (in deg) for track candidate clustering |
312 | // into jets. | |
67338e65 | 313 | // In the constructor the default has been set to 15 deg, in accordance |
314 | // to NIM A524 (2004) 180. | |
315 | // | |
413d0114 | 316 | // Note : This function also sets automatically the value of the maximum |
317 | // angular separation for jet merging into 1 single track to ang/2. | |
318 | // In order to specify a different max. jet merging separation angle, | |
319 | // one has to invoke the memberfunction SetJangmax afterwards. | |
320 | ||
321 | fTangmax=ang; | |
322 | fJangmax=ang/2.; | |
323 | } | |
324 | /////////////////////////////////////////////////////////////////////////// | |
25eefd00 | 325 | void IceDwalk::SetTdistmax(Float_t d,Int_t invol) |
413d0114 | 326 | { |
25eefd00 | 327 | // Set maximum distance (in m) of the two track candidates in the track |
328 | // clustering process. | |
329 | // The distance between the two tracks can be determined restricted to the | |
330 | // detector volume (invol=1) or in the overall space (invol=0). | |
331 | // The former will prevent clustering of (nearly) parallel tracks which cross | |
332 | // the detector volume at very different locations, whereas the latter will | |
333 | // enable clustering of tracks with a common location of origin (e.g. muon | |
334 | // bundles from an air shower) even if they cross the detector volume at | |
335 | // very different locations. | |
336 | // At invokation of this memberfunction the default is invol=1. | |
337 | // In the constructor the default has been set to 20 meter with invol=1. | |
67338e65 | 338 | |
25eefd00 | 339 | fTdistmax=d; |
340 | fTinvol=invol; | |
67338e65 | 341 | } |
342 | /////////////////////////////////////////////////////////////////////////// | |
25eefd00 | 343 | void IceDwalk::SetJangmax(Float_t ang,Int_t iter) |
67338e65 | 344 | { |
345 | // Set angular separation (in deg) within which jets are merged into 1 | |
346 | // single track. | |
25eefd00 | 347 | // The merging process is a dynamic procedure and can be carried out by |
348 | // iteration (iter=1) until no further merging of the various jets occurs anymore. | |
349 | // However, by specification of iter=0 the user can also select to go only | |
350 | // once through all the jet combinations to check for mergers. | |
351 | // For large events the latter will in general result in more track candidates. | |
352 | // At invokation of this memberfunction the default is iter=1. | |
353 | // In the constructor the default angle has been set 7.5 deg, being half | |
354 | // of the value of the default track candidate clustering separation angle. | |
355 | // The iteration flag was set to 1 in the constructor. | |
67338e65 | 356 | // |
357 | // Notes : | |
358 | // ------- | |
359 | // 1) Setting ang=0 will prevent jet merging. | |
360 | // Consequently, every jet will appear as a separate track in the | |
361 | // reconstruction result. | |
362 | // 2) Setting ang<0 will prevent jet merging. | |
363 | // In addition, only the jet with the maximum number of tracks will | |
364 | // appear as a track in the reconstruction result. | |
365 | // This situation resembles the standard Sieglinde direct walk processing | |
366 | // and as such can be used to perform comparison studies. | |
367 | ||
413d0114 | 368 | fJangmax=ang; |
25eefd00 | 369 | fJiterate=iter; |
413d0114 | 370 | } |
371 | /////////////////////////////////////////////////////////////////////////// | |
25eefd00 | 372 | void IceDwalk::SetJdistmax(Float_t d,Int_t invol) |
413d0114 | 373 | { |
25eefd00 | 374 | // Set maximum distance (in m) of the two jets in the jet merging process. |
375 | // The distance between the two jets can be determined restricted to the | |
376 | // detector volume (invol=1) or in the overall space (invol=0). | |
377 | // The former will prevent clustering of (nearly) parallel tracks which cross | |
378 | // the detector volume at very different locations, whereas the latter will | |
379 | // enable clustering of tracks with a common location of origin (e.g. muon | |
380 | // bundles from an air shower) even if they cross the detector volume at | |
381 | // very different locations. | |
382 | // At invokation of this memberfunction the default is invol=1. | |
383 | // In the constructor the default has been set to 30 meter with invol=1. | |
413d0114 | 384 | |
25eefd00 | 385 | fJdistmax=d; |
386 | fJinvol=invol; | |
67338e65 | 387 | } |
388 | /////////////////////////////////////////////////////////////////////////// | |
d6860cf1 | 389 | void IceDwalk::SetMaxModA(Int_t nmax) |
390 | { | |
391 | // Set the maximum number of good Amanda modules that may have fired | |
392 | // in order to process this event. | |
393 | // This allows suppression of processing (high-energy) cascade events | |
394 | // with this direct walk tracking to prevent waisting cpu time for cases | |
395 | // in which tracking doesn't make sense anyhow. | |
396 | // Furthermore it allows selection of low multiplicity events for processing. | |
397 | // By default the maximum number of Amanda modules is set to 999 in the ctor, | |
6f94f699 | 398 | // which implies no selection on maximum module multiplicity. |
399 | // See also the memberfunction SetMinModA(). | |
d6860cf1 | 400 | fMaxmodA=nmax; |
401 | } | |
402 | /////////////////////////////////////////////////////////////////////////// | |
6f94f699 | 403 | void IceDwalk::SetMinModA(Int_t nmin) |
404 | { | |
405 | // Set the minimum number of good Amanda modules that must have fired | |
406 | // in order to process this event. | |
407 | // This allows selection of a minimal multiplicity for events to be processed. | |
408 | // By default the minimum number of Amanda modules is set to 0 in the ctor, | |
409 | // which implies no selection on minimum module multiplicity. | |
410 | // See also the memberfunction SetMaxModA(). | |
411 | fMinmodA=nmin; | |
412 | } | |
413 | /////////////////////////////////////////////////////////////////////////// | |
35e721bf | 414 | void IceDwalk::SetMaxHitsA(Int_t nmax) |
415 | { | |
25eefd00 | 416 | // Set the maximum number of good hits per Amanda module to be processed. |
35e721bf | 417 | // |
418 | // Special values : | |
25eefd00 | 419 | // nmax = 0 : No maximum limit set; all good hits will be processed |
35e721bf | 420 | // < 0 : No hits will be processed |
421 | // | |
25eefd00 | 422 | // In case the user selects a maximum number of good hits per module, all the |
35e721bf | 423 | // hits of each module will be ordered w.r.t. increasing hit time (LE). |
424 | // This allows selection of processing e.g. only the first hits etc... | |
25eefd00 | 425 | // By default the maximum number of hits per Amanda modules is set to 1 in the ctor, |
426 | // which implies processing only the first good hit of each Amanda OM. | |
35e721bf | 427 | fMaxhitsA=nmax; |
428 | } | |
429 | /////////////////////////////////////////////////////////////////////////// | |
25eefd00 | 430 | void IceDwalk::SetVgroupUsage(Int_t flag) |
431 | { | |
432 | // (De)activate the distinction between v_phase and v_group of the Cherenkov light. | |
433 | // | |
434 | // flag = 0 : No distinction between v_phase and v_group | |
435 | // = 1 : Separate treatment of v_phase and v_group | |
436 | // | |
437 | // By default the distinction between v_phase and v_group is activated | |
438 | // in the constructor of this class. | |
439 | fVgroup=flag; | |
440 | } | |
441 | /////////////////////////////////////////////////////////////////////////// | |
d6860cf1 | 442 | void IceDwalk::SetTrackName(TString s) |
443 | { | |
444 | // Set (alternative) name identifier for the produced first guess tracks. | |
445 | // This allows unique identification of (newly) produced direct walk tracks | |
446 | // in case of re-processing of existing data with different criteria. | |
bfde1ea0 | 447 | // By default the produced first guess tracks have the name of the class |
448 | // by which they were produced. | |
d6860cf1 | 449 | fTrackname=s; |
450 | } | |
451 | /////////////////////////////////////////////////////////////////////////// | |
30672a96 | 452 | void IceDwalk::SetCharge(Float_t charge) |
453 | { | |
454 | // Set user defined charge for the produced first guess tracks. | |
455 | // This allows identification of these tracks on color displays. | |
456 | // By default the produced first guess tracks have charge=0 | |
457 | // which is set in the constructor of this class. | |
458 | fCharge=charge; | |
459 | } | |
460 | /////////////////////////////////////////////////////////////////////////// | |
67338e65 | 461 | void IceDwalk::Exec(Option_t* opt) |
462 | { | |
463 | // Implementation of the direct walk track reconstruction. | |
464 | ||
465 | TString name=opt; | |
466 | AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data())); | |
467 | ||
468 | if (!parent) return; | |
469 | ||
bfde1ea0 | 470 | fEvt=(IceEvent*)parent->GetObject("IceEvent"); |
471 | if (!fEvt) return; | |
67338e65 | 472 | |
caa58e1a | 473 | // Only process accepted events |
474 | AliDevice* seldev=(AliDevice*)fEvt->GetDevice("AliEventSelector"); | |
475 | if (seldev) | |
476 | { | |
477 | if (seldev->GetSignal("Select") < 0.1) return; | |
478 | } | |
479 | ||
35e721bf | 480 | // Enter the reco parameters as a device in the event |
481 | AliSignal params; | |
bfde1ea0 | 482 | params.SetNameTitle(ClassName(),"Reco parameters"); |
35e721bf | 483 | params.SetSlotName("Dmin",1); |
484 | params.SetSlotName("Dtmarg",2); | |
25eefd00 | 485 | params.SetSlotName("Maxdhit",3); |
486 | params.SetSlotName("Tangmax",4); | |
487 | params.SetSlotName("Tdistmax",5); | |
488 | params.SetSlotName("Tinvol",6); | |
489 | params.SetSlotName("Jangmax",7); | |
490 | params.SetSlotName("Jiterate",8); | |
491 | params.SetSlotName("Jdistmax",9); | |
492 | params.SetSlotName("Jinvol",10); | |
493 | params.SetSlotName("MaxmodA",11); | |
494 | params.SetSlotName("MinmodA",12); | |
495 | params.SetSlotName("MaxhitsA",13); | |
496 | params.SetSlotName("Vgroup",14); | |
35e721bf | 497 | |
498 | params.SetSignal(fDmin,1); | |
499 | params.SetSignal(fDtmarg,2); | |
25eefd00 | 500 | params.SetSignal(fMaxdhit,3); |
501 | params.SetSignal(fTangmax,4); | |
502 | params.SetSignal(fTdistmax,5); | |
503 | params.SetSignal(fTinvol,6); | |
504 | params.SetSignal(fJangmax,7); | |
505 | params.SetSignal(fJiterate,8); | |
506 | params.SetSignal(fJdistmax,9); | |
507 | params.SetSignal(fJinvol,10); | |
508 | params.SetSignal(fMaxmodA,11); | |
509 | params.SetSignal(fMinmodA,12); | |
510 | params.SetSignal(fMaxhitsA,13); | |
511 | params.SetSignal(fVgroup,14); | |
35e721bf | 512 | |
bfde1ea0 | 513 | fEvt->AddDevice(params); |
35e721bf | 514 | |
515 | if (fMaxhitsA<0) return; | |
516 | ||
413d0114 | 517 | // Fetch all fired Amanda OMs for this event |
bfde1ea0 | 518 | TObjArray* aoms=fEvt->GetDevices("IceAOM"); |
8a394508 | 519 | if (!aoms) return; |
413d0114 | 520 | Int_t naoms=aoms->GetEntries(); |
521 | if (!naoms) return; | |
522 | ||
6f94f699 | 523 | // Check for the minimum and/or maximum number of good fired Amanda OMs |
d6860cf1 | 524 | Int_t ngood=0; |
525 | for (Int_t iom=0; iom<naoms; iom++) | |
526 | { | |
527 | IceGOM* omx=(IceGOM*)aoms->At(iom); | |
528 | if (!omx) continue; | |
529 | if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue; | |
530 | ngood++; | |
531 | } | |
6f94f699 | 532 | if (ngood<fMinmodA || ngood>fMaxmodA) return; |
d6860cf1 | 533 | |
bfde1ea0 | 534 | const Float_t c=0.299792458; // Light speed in vacuum in meters per ns |
67338e65 | 535 | |
35e721bf | 536 | // Storage of track elements. |
67338e65 | 537 | TObjArray tes; |
538 | tes.SetOwner(); | |
67338e65 | 539 | |
540 | AliPosition r1; | |
541 | AliPosition r2; | |
542 | Ali3Vector r12; | |
80701e9a | 543 | Ali3Vector rsum; |
67338e65 | 544 | AliPosition r0; |
545 | TObjArray hits1; | |
546 | TObjArray hits2; | |
547 | Int_t nh1,nh2; | |
548 | AliSignal* sx1=0; | |
549 | AliSignal* sx2=0; | |
c8d547c3 | 550 | Float_t dist=0; |
67338e65 | 551 | Float_t t1,t2,dt,t0; |
35e721bf | 552 | Float_t dtmax; |
67338e65 | 553 | TObjArray hits; |
35e721bf | 554 | TObjArray* ordered; |
67338e65 | 555 | |
bfde1ea0 | 556 | // Check the hits of Amanda OM pairs for possible track elements. |
67338e65 | 557 | // Also all the good hits are stored in the meantime (to save CPU time) |
558 | // for hit association with the various track elements lateron. | |
67338e65 | 559 | AliTrack* te=0; |
67338e65 | 560 | for (Int_t i1=0; i1<naoms; i1++) // First OM of the pair |
561 | { | |
562 | IceGOM* omx1=(IceGOM*)aoms->At(i1); | |
563 | if (!omx1) continue; | |
564 | if (omx1->GetDeadValue("LE")) continue; | |
565 | r1=omx1->GetPosition(); | |
566 | // Select all the good hits of this first OM | |
567 | hits1.Clear(); | |
35e721bf | 568 | // Determine the max. number of hits to be processed for this OM |
569 | ordered=0; | |
570 | if (fMaxhitsA>0 && omx1->GetNhits()>fMaxhitsA) ordered=omx1->SortHits("LE",1,0,7); | |
571 | nh1=0; | |
67338e65 | 572 | for (Int_t j1=1; j1<=omx1->GetNhits(); j1++) |
573 | { | |
35e721bf | 574 | if (ordered) |
575 | { | |
576 | if (nh1>=fMaxhitsA) break; | |
577 | sx1=(AliSignal*)ordered->At(j1-1); | |
578 | } | |
579 | else | |
580 | { | |
581 | sx1=omx1->GetHit(j1); | |
582 | } | |
67338e65 | 583 | if (!sx1) continue; |
584 | if (sx1->GetDeadValue("ADC") || sx1->GetDeadValue("LE") || sx1->GetDeadValue("TOT")) continue; | |
585 | hits1.Add(sx1); | |
586 | // Also store all good hits in the total hit array | |
587 | hits.Add(sx1); | |
35e721bf | 588 | nh1++; |
67338e65 | 589 | } |
590 | ||
591 | // No further pair to be formed with the last OM in the list | |
592 | if (i1==(naoms-1)) break; | |
593 | ||
594 | nh1=hits1.GetEntries(); | |
595 | if (!nh1) continue; | |
596 | ||
597 | for (Int_t i2=i1+1; i2<naoms; i2++) // Second OM of the pair | |
598 | { | |
599 | IceGOM* omx2=(IceGOM*)aoms->At(i2); | |
600 | if (!omx2) continue; | |
601 | if (omx2->GetDeadValue("LE")) continue; | |
602 | r2=omx2->GetPosition(); | |
603 | r12=r2-r1; | |
604 | dist=r12.GetNorm(); | |
413d0114 | 605 | |
67338e65 | 606 | if (dist<=fDmin) continue; |
413d0114 | 607 | |
67338e65 | 608 | // Select all the good hits of this second OM |
609 | hits2.Clear(); | |
35e721bf | 610 | // Determine the max. number of hits to be processed for this OM |
611 | ordered=0; | |
612 | if (fMaxhitsA>0 && omx2->GetNhits()>fMaxhitsA) ordered=omx2->SortHits("LE",1,0,7); | |
613 | nh2=0; | |
67338e65 | 614 | for (Int_t j2=1; j2<=omx2->GetNhits(); j2++) |
615 | { | |
35e721bf | 616 | if (ordered) |
617 | { | |
618 | if (nh2>=fMaxhitsA) break; | |
619 | sx2=(AliSignal*)ordered->At(j2-1); | |
620 | } | |
621 | else | |
622 | { | |
623 | sx2=omx2->GetHit(j2); | |
624 | } | |
67338e65 | 625 | if (!sx2) continue; |
626 | if (sx2->GetDeadValue("ADC") || sx2->GetDeadValue("LE") || sx2->GetDeadValue("TOT")) continue; | |
627 | hits2.Add(sx2); | |
35e721bf | 628 | nh2++; |
67338e65 | 629 | } |
630 | ||
631 | nh2=hits2.GetEntries(); | |
632 | if (!nh2) continue; | |
633 | ||
413d0114 | 634 | // Position r0 in between the two OMs and normalised relative direction r12 |
635 | rsum=(r1+r2)/2.; | |
636 | r0.SetPosition((Ali3Vector&)rsum); | |
637 | r12/=dist; | |
638 | ||
67338e65 | 639 | // Check all hit pair combinations of these two OMs for possible track elements |
413d0114 | 640 | dtmax=dist/c+float(fDtmarg); |
67338e65 | 641 | for (Int_t ih1=0; ih1<nh1; ih1++) // Hits of first OM |
642 | { | |
643 | sx1=(AliSignal*)hits1.At(ih1); | |
644 | if (!sx1) continue; | |
645 | for (Int_t ih2=0; ih2<nh2; ih2++) // Hits of second OM | |
646 | { | |
647 | sx2=(AliSignal*)hits2.At(ih2); | |
648 | if (!sx2) continue; | |
649 | t1=sx1->GetSignal("LE",7); | |
650 | t2=sx2->GetSignal("LE",7); | |
651 | dt=t2-t1; | |
652 | t0=(t1+t2)/2.; | |
413d0114 | 653 | |
654 | if (fabs(dt)>=dtmax) continue; | |
655 | ||
656 | te=new AliTrack(); | |
657 | tes.Add(te); | |
413d0114 | 658 | if (dt<0) r12*=-1.; |
bfde1ea0 | 659 | r0.SetTimestamp((AliTimestamp&)*fEvt); |
413d0114 | 660 | AliTimestamp* tsx=r0.GetTimestamp(); |
661 | tsx->Add(0,0,(int)t0); | |
662 | te->SetReferencePoint(r0); | |
663 | te->Set3Momentum(r12); | |
67338e65 | 664 | } |
665 | } | |
666 | } // end of loop over the second OM of the pair | |
667 | } // end of loop over first OM of the pair | |
668 | ||
25eefd00 | 669 | // Association of hits to the various track elements. |
bfde1ea0 | 670 | Float_t qmax=0; |
671 | Int_t nahmax=0; | |
672 | AssociateHits(tes,hits,qmax,nahmax); | |
673 | ||
674 | // Selection on quality (Q value) in case of multiple track candidates | |
675 | SelectQvalue(tes,qmax); | |
676 | ||
677 | Int_t nte=tes.GetEntries(); | |
678 | if (!nte) return; | |
679 | ||
680 | // Clustering of track candidates into jets | |
681 | TObjArray jets; | |
682 | jets.SetOwner(); | |
683 | ClusterTracks(tes,jets,qmax); | |
684 | ||
685 | Int_t njets=jets.GetEntries(); | |
686 | if (!njets) return; | |
687 | ||
688 | // Order the jets w.r.t. decreasing quality value | |
689 | ordered=fEvt->SortJets(-2,&jets); | |
690 | TObjArray jets2(*ordered); | |
691 | ||
692 | // Merging f jets | |
693 | MergeJets(jets2); | |
694 | ||
695 | // Production and storage of the final tracks | |
696 | StoreTracks(jets2); | |
697 | } | |
698 | /////////////////////////////////////////////////////////////////////////// | |
699 | void IceDwalk::AssociateHits(TObjArray& tes,TObjArray& hits,Float_t& qmax,Int_t& nahmax) | |
700 | { | |
701 | // Association of hits to the various track elements. | |
702 | ||
703 | const Float_t pi=acos(-1.); | |
704 | const Float_t c=0.299792458; // Light speed in vacuum in meters per ns | |
705 | const Float_t npice=1.31768387; // Phase refractive index (c/v_phase) of ice | |
706 | const Float_t ngice=1.35075806; // Group refractive index (c/v_group) of ice | |
707 | const Float_t lambda=33.3; // Light scattering length in ice | |
708 | const Float_t thetac=acos(1./npice); // Cherenkov angle (in radians) | |
709 | ||
710 | // Angular reduction of complement of thetac due to v_phase and v_group difference | |
711 | Float_t alphac=0; | |
712 | if (fVgroup) alphac=atan((1.-npice/ngice)/sqrt(npice*npice-1.)); | |
713 | ||
67338e65 | 714 | Int_t nte=tes.GetEntries(); |
715 | Int_t nh=hits.GetEntries(); | |
716 | Float_t d=0; | |
717 | Ali3Vector p; | |
bfde1ea0 | 718 | AliPosition r1; |
719 | AliPosition r2; | |
720 | Ali3Vector r12; | |
721 | Float_t t1; | |
722 | Float_t dist,t0,tgeo,tres; | |
25eefd00 | 723 | AliSample levers; // Statistics of the assoc. hit lever arms |
724 | levers.SetStoreMode(1);// Enable median calculation | |
725 | AliSample hprojs; // Statistics of the assoc. hit position projections on the track w.r.t. r0 | |
726 | hprojs.SetStoreMode(1);// Enable median calculation | |
727 | AliSample times; // Statistics of the time residuals of the associated hits | |
728 | times.SetStoreMode(1); // Enable median calculation | |
729 | AliSignal fit; // Storage of Q value etc... for each track candidate | |
730 | fit.AddNamedSlot("QTC"); | |
731 | fit.AddNamedSlot("SpanL"); | |
732 | fit.AddNamedSlot("MedianL"); | |
733 | fit.AddNamedSlot("MeanL"); | |
734 | fit.AddNamedSlot("SigmaL"); | |
735 | fit.AddNamedSlot("SpreadL"); | |
736 | fit.AddNamedSlot("ExpSpreadL"); | |
737 | fit.AddNamedSlot("Span"); | |
738 | fit.AddNamedSlot("Median"); | |
739 | fit.AddNamedSlot("Mean"); | |
740 | fit.AddNamedSlot("Sigma"); | |
741 | fit.AddNamedSlot("Spread"); | |
742 | fit.AddNamedSlot("ExpSpread"); | |
743 | fit.AddNamedSlot("MedianT"); | |
744 | fit.AddNamedSlot("MeanT"); | |
745 | fit.AddNamedSlot("SigmaT"); | |
746 | fit.AddNamedSlot("SpreadT"); | |
747 | fit.AddNamedSlot("term1"); | |
748 | fit.AddNamedSlot("term2"); | |
749 | fit.AddNamedSlot("term3"); | |
750 | fit.AddNamedSlot("term4"); | |
751 | fit.AddNamedSlot("term5"); | |
bfde1ea0 | 752 | Float_t qtc=0; |
753 | Int_t nah; // Number of associated hits for a certain TE | |
25eefd00 | 754 | Float_t lmin,lmax,spanl,medianl,meanl,sigmal,spreadl,expspreadl; |
755 | Float_t hproj,hprojmin,hprojmax,span,median,mean,sigma,spread,expspread; | |
756 | Float_t mediant,meant,sigmat,spreadt; | |
757 | Float_t term1,term2,term3,term4,term5; | |
bfde1ea0 | 758 | qmax=0; |
759 | nahmax=0; | |
67338e65 | 760 | for (Int_t jte=0; jte<nte; jte++) |
761 | { | |
bfde1ea0 | 762 | AliTrack* te=(AliTrack*)tes.At(jte); |
67338e65 | 763 | if (!te) continue; |
67338e65 | 764 | AliPosition* tr0=te->GetReferencePoint(); |
413d0114 | 765 | AliTimestamp* tt0=tr0->GetTimestamp(); |
bfde1ea0 | 766 | t0=fEvt->GetDifference(tt0,"ns"); |
413d0114 | 767 | p=te->Get3Momentum(); |
67338e65 | 768 | levers.Reset(); |
25eefd00 | 769 | hprojs.Reset(); |
770 | times.Reset(); | |
67338e65 | 771 | for (Int_t jh=0; jh<nh; jh++) |
772 | { | |
bfde1ea0 | 773 | AliSignal* sx1=(AliSignal*)hits.At(jh); |
67338e65 | 774 | if (!sx1) continue; |
775 | IceGOM* omx=(IceGOM*)sx1->GetDevice(); | |
776 | if (!omx) continue; | |
777 | r1=omx->GetPosition(); | |
35e721bf | 778 | d=te->GetDistance(r1); |
67338e65 | 779 | r12=r1-(*tr0); |
25eefd00 | 780 | hproj=p.Dot(r12); |
781 | dist=hproj+d/tan(pi/2.-thetac-alphac); | |
67338e65 | 782 | tgeo=t0+dist/c; |
783 | t1=sx1->GetSignal("LE",7); | |
784 | tres=t1-tgeo; | |
785 | ||
bfde1ea0 | 786 | d=d/sin(thetac); // The distance traveled by a cherenkov photon |
787 | ||
25eefd00 | 788 | if (tres<-30 || tres>300 || d>fMaxdhit*lambda) continue; |
67338e65 | 789 | |
790 | // Associate this hit to the TE | |
791 | te->AddSignal(*sx1); | |
25eefd00 | 792 | levers.Enter(fabs(hproj)); |
793 | hprojs.Enter(hproj); | |
794 | times.Enter(tres); | |
67338e65 | 795 | } |
413d0114 | 796 | |
797 | // Determine the Q quality of the various TE's. | |
798 | // Good quality TE's will be called track candidates (TC's) | |
67338e65 | 799 | nah=te->GetNsignals(); |
25eefd00 | 800 | if (nah>nahmax) nahmax=nah; |
801 | lmin=levers.GetMinimum(1); | |
802 | lmax=levers.GetMaximum(1); | |
803 | spanl=lmax-lmin; | |
804 | medianl=levers.GetMedian(1); | |
805 | meanl=levers.GetMean(1); | |
67338e65 | 806 | sigmal=levers.GetSigma(1); |
25eefd00 | 807 | spreadl=levers.GetSpread(1); |
808 | // Expected spread for a flat distribution | |
809 | expspreadl=0; | |
810 | if (spanl>0) expspreadl=(0.5*pow(lmin,2)+0.5*pow(lmax,2)+pow(medianl,2)-medianl*(lmin+lmax))/spanl; | |
811 | hprojmin=hprojs.GetMinimum(1); | |
812 | hprojmax=hprojs.GetMaximum(1); | |
813 | span=hprojmax-hprojmin; | |
814 | median=hprojs.GetMedian(1); | |
815 | mean=hprojs.GetMean(1); | |
816 | sigma=hprojs.GetSigma(1); | |
817 | spread=hprojs.GetSpread(1); | |
818 | // Expected spread for a flat distribution | |
819 | expspread=0; | |
820 | if (span>0) expspread=(0.5*pow(hprojmin,2)+0.5*pow(hprojmax,2)+pow(median,2)-median*(hprojmin+hprojmax))/span; | |
821 | mediant=times.GetMedian(1); | |
822 | meant=times.GetMean(1); | |
823 | sigmat=times.GetSigma(1); | |
824 | spreadt=times.GetSpread(1); | |
825 | ||
826 | term1=0; | |
827 | if (span>0) term1=2.*spread/span; | |
828 | ||
829 | term2=0; | |
830 | if (spanl>0) term2=2.*spreadl/spanl; | |
831 | ||
832 | term3=0; | |
833 | if (spread>0) term3=fabs(spread-expspread)/spread; | |
834 | ||
835 | term4=0; | |
836 | if (spreadl>0) term4=fabs(spreadl-expspreadl)/spreadl; | |
837 | ||
838 | term5=0; | |
839 | if (spreadt>0) term5=fabs(mediant)/spreadt; | |
840 | ||
841 | qtc=float(nah)*(term1+term2)-term3-term4-term5; | |
842 | if (fabs(median)>span/2.) qtc=0; // Require projected hits on both sides of r0 | |
843 | ||
413d0114 | 844 | fit.SetSignal(qtc,"QTC"); |
25eefd00 | 845 | fit.SetSignal(spanl,"SpanL"); |
846 | fit.SetSignal(medianl,"MedianL"); | |
847 | fit.SetSignal(meanl,"MeanL"); | |
848 | fit.SetSignal(sigmal,"SigmaL"); | |
849 | fit.SetSignal(spreadl,"SpreadL"); | |
850 | fit.SetSignal(expspreadl,"ExpSpreadL"); | |
851 | fit.SetSignal(span,"Span"); | |
852 | fit.SetSignal(median,"Median"); | |
853 | fit.SetSignal(mean,"Mean"); | |
854 | fit.SetSignal(sigma,"Sigma"); | |
855 | fit.SetSignal(spread,"Spread"); | |
856 | fit.SetSignal(expspread,"ExpSpread"); | |
857 | fit.SetSignal(mediant,"MedianT"); | |
858 | fit.SetSignal(meant,"MeanT"); | |
859 | fit.SetSignal(sigmat,"SigmaT"); | |
860 | fit.SetSignal(spreadt,"SpreadT"); | |
861 | fit.SetSignal(term1,"term1"); | |
862 | fit.SetSignal(term2,"term2"); | |
863 | fit.SetSignal(term3,"term3"); | |
864 | fit.SetSignal(term4,"term4"); | |
865 | fit.SetSignal(term5,"term5"); | |
413d0114 | 866 | te->SetFitDetails(fit); |
867 | if (qtc>qmax) qmax=qtc; | |
67338e65 | 868 | } |
bfde1ea0 | 869 | } |
870 | /////////////////////////////////////////////////////////////////////////// | |
871 | void IceDwalk::SelectQvalue(TObjArray& tes,Float_t qmax) | |
872 | { | |
67338e65 | 873 | // Perform selection on Q value in case of multiple track candidates |
bfde1ea0 | 874 | |
875 | Int_t nte=tes.GetEntries(); | |
876 | Int_t nah; | |
877 | Float_t qtc; | |
878 | Ali3Vector p; | |
67338e65 | 879 | for (Int_t jtc=0; jtc<nte; jtc++) |
880 | { | |
bfde1ea0 | 881 | AliTrack* te=(AliTrack*)tes.At(jtc); |
67338e65 | 882 | if (!te) continue; |
25eefd00 | 883 | nah=te->GetNsignals(); |
bfde1ea0 | 884 | AliSignal* sx1=(AliSignal*)te->GetFitDetails(); |
413d0114 | 885 | qtc=-1; |
886 | if (sx1) qtc=sx1->GetSignal("QTC"); | |
25eefd00 | 887 | if (!nah || qtc<0.8*qmax) |
67338e65 | 888 | { |
67338e65 | 889 | tes.RemoveAt(jtc); |
890 | delete te; | |
891 | } | |
25eefd00 | 892 | else // Set Q value as momentum to provide a weight for jet clustering |
893 | { | |
894 | if (qtc>0) | |
895 | { | |
896 | p=te->Get3Momentum(); | |
897 | p*=qtc; | |
898 | te->Set3Momentum(p); | |
899 | } | |
900 | } | |
67338e65 | 901 | } |
902 | tes.Compress(); | |
bfde1ea0 | 903 | } |
904 | /////////////////////////////////////////////////////////////////////////// | |
905 | void IceDwalk::ClusterTracks(TObjArray& tes,TObjArray& jets,Float_t qmax) | |
906 | { | |
413d0114 | 907 | // Cluster track candidates within a certain opening angle into jets. |
25eefd00 | 908 | // Also the track should be within a certain maximum distance of the |
909 | // starting track in order to get clustered. | |
413d0114 | 910 | // The latter prevents clustering of (nearly) parallel track candidates |
911 | // crossing the detector a very different locations (e.g. muon bundles). | |
912 | // The average r0 and t0 of the constituent tracks will be taken as the | |
913 | // jet reference point. | |
bfde1ea0 | 914 | |
915 | Int_t nte=tes.GetEntries(); | |
67338e65 | 916 | Float_t ang=0; |
413d0114 | 917 | AliSample pos; |
918 | AliSample time; | |
919 | Float_t vec[3],err[3]; | |
bfde1ea0 | 920 | AliPosition r0; |
921 | Float_t t0,dist,dist2; | |
922 | Int_t nah=0,nahmax=0; // Determine the max. number of associated hits for the jets | |
923 | Float_t qtc; | |
67338e65 | 924 | for (Int_t jtc1=0; jtc1<nte; jtc1++) |
925 | { | |
bfde1ea0 | 926 | AliTrack* te=(AliTrack*)tes.At(jtc1); |
67338e65 | 927 | if (!te) continue; |
413d0114 | 928 | AliPosition* x1=te->GetReferencePoint(); |
929 | if (!x1) continue; | |
930 | AliTimestamp* ts1=x1->GetTimestamp(); | |
931 | if (!ts1) continue; | |
67338e65 | 932 | AliJet* jx=new AliJet(); |
933 | jx->AddTrack(te); | |
413d0114 | 934 | pos.Reset(); |
935 | time.Reset(); | |
936 | x1->GetPosition(vec,"car"); | |
937 | pos.Enter(vec[0],vec[1],vec[2]); | |
bfde1ea0 | 938 | t0=fEvt->GetDifference(ts1,"ns"); |
413d0114 | 939 | time.Enter(t0); |
67338e65 | 940 | for (Int_t jtc2=0; jtc2<nte; jtc2++) |
941 | { | |
413d0114 | 942 | if (jtc2==jtc1) continue; |
bfde1ea0 | 943 | AliTrack* te2=(AliTrack*)tes.At(jtc2); |
67338e65 | 944 | if (!te2) continue; |
945 | ang=te->GetOpeningAngle(*te2,"deg"); | |
25eefd00 | 946 | if (ang<=fTangmax) |
413d0114 | 947 | { |
948 | AliPosition* x2=te2->GetReferencePoint(); | |
949 | if (!x2) continue; | |
950 | AliTimestamp* ts2=x2->GetTimestamp(); | |
951 | if (!ts2) continue; | |
25eefd00 | 952 | if (!fTinvol) |
413d0114 | 953 | { |
25eefd00 | 954 | dist=te->GetDistance(te2); |
955 | } | |
956 | else | |
957 | { | |
958 | dist=te->GetDistance(x2); | |
959 | dist2=te2->GetDistance(x1); | |
960 | if (dist2<dist) dist=dist2; | |
961 | } | |
962 | if (dist<=fTdistmax) | |
963 | { | |
964 | x2->GetPosition(vec,"car"); | |
965 | pos.Enter(vec[0],vec[1],vec[2]); | |
bfde1ea0 | 966 | t0=fEvt->GetDifference(ts2,"ns"); |
25eefd00 | 967 | time.Enter(t0); |
968 | jx->AddTrack(te2); | |
413d0114 | 969 | } |
970 | } | |
971 | } | |
972 | ||
973 | // Set the reference point data for this jet | |
974 | for (Int_t j=1; j<=3; j++) | |
975 | { | |
976 | vec[j-1]=pos.GetMean(j); | |
977 | err[j-1]=pos.GetSigma(j); | |
978 | } | |
979 | r0.SetPosition(vec,"car"); | |
980 | r0.SetPositionErrors(err,"car"); | |
bfde1ea0 | 981 | r0.SetTimestamp((AliTimestamp&)*fEvt); |
413d0114 | 982 | AliTimestamp* jt0=r0.GetTimestamp(); |
983 | t0=time.GetMean(1); | |
984 | jt0->Add(0,0,(int)t0); | |
985 | jx->SetReferencePoint(r0); | |
986 | ||
987 | // Store this jet for further processing if ntracks>1 | |
25eefd00 | 988 | if (jx->GetNtracks() > 1 || fTangmax<=0) |
413d0114 | 989 | { |
990 | jets.Add(jx); | |
25eefd00 | 991 | nah=jx->GetNsignals(); |
992 | if (nah>nahmax) nahmax=nah; | |
413d0114 | 993 | } |
994 | else // Only keep single-track jets which have qtc=qmax | |
995 | { | |
bfde1ea0 | 996 | AliSignal* sx1=(AliSignal*)te->GetFitDetails(); |
413d0114 | 997 | qtc=-1; |
998 | if (sx1) qtc=sx1->GetSignal("QTC"); | |
999 | if (qtc>=(qmax-1.e-10)) | |
1000 | { | |
1001 | jets.Add(jx); | |
25eefd00 | 1002 | nah=jx->GetNsignals(); |
1003 | if (nah>nahmax) nahmax=nah; | |
413d0114 | 1004 | } |
1005 | else | |
1006 | { | |
1007 | delete jx; | |
1008 | } | |
67338e65 | 1009 | } |
1010 | } | |
413d0114 | 1011 | |
bfde1ea0 | 1012 | Int_t njets=jets.GetEntries(); |
413d0114 | 1013 | if (!njets) return; |
67338e65 | 1014 | |
bfde1ea0 | 1015 | // The sum of 0.15*(nah-nahmax) and average qtc value per track for each jet |
25eefd00 | 1016 | // will be stored as the jet energy to enable sorting on this value lateron |
1017 | Float_t sortval=0; | |
1018 | Int_t ntk=0; | |
1019 | for (Int_t ijet=0; ijet<njets; ijet++) | |
1020 | { | |
1021 | AliJet* jx=(AliJet*)jets.At(ijet); | |
1022 | if (!jx) continue; | |
1023 | nah=jx->GetNsignals(); | |
1024 | ntk=jx->GetNtracks(); | |
1025 | sortval=0.15*float(nah-nahmax); | |
1026 | if (ntk) sortval+=jx->GetMomentum()/float(ntk); | |
1027 | jx->SetScalar(sortval); | |
1028 | } | |
bfde1ea0 | 1029 | } |
1030 | /////////////////////////////////////////////////////////////////////////// | |
1031 | void IceDwalk::MergeJets(TObjArray& jets2) | |
1032 | { | |
67338e65 | 1033 | // Merge jets within a certain opening to provide the final track(s). |
25eefd00 | 1034 | // Also the jet should be within a certain maximum distance of the |
1035 | // starting jet in order to get merged. | |
1036 | // The latter prevents merging of (nearly) parallel jets/tracks | |
1037 | // crossing the detector a very different locations (e.g. muon bundles). | |
1038 | // The average r0 and t0 of the constituent jets will be taken as the | |
1039 | // final reference point. | |
bfde1ea0 | 1040 | |
1041 | Int_t njets=jets2.GetEntries(); | |
67338e65 | 1042 | AliJet* jx1=0; |
1043 | AliJet* jx2=0; | |
25eefd00 | 1044 | Int_t merged=1; |
bfde1ea0 | 1045 | Int_t ntk,nah,nahmax; |
1046 | Float_t ang,dist,dist2,t0; | |
1047 | AliSample pos; | |
1048 | AliSample time; | |
1049 | AliPosition r0; | |
1050 | Float_t vec[3],err[3]; | |
1051 | Float_t sortval; | |
25eefd00 | 1052 | if (fJangmax>=0) |
67338e65 | 1053 | { |
25eefd00 | 1054 | while (merged) |
67338e65 | 1055 | { |
25eefd00 | 1056 | merged=0; |
1057 | nahmax=0; | |
1058 | for (Int_t jet1=0; jet1<njets; jet1++) | |
67338e65 | 1059 | { |
25eefd00 | 1060 | jx1=(AliJet*)jets2.At(jet1); |
1061 | if (!jx1) continue; | |
1062 | AliPosition* x1=jx1->GetReferencePoint(); | |
1063 | if (!x1) continue; | |
1064 | AliTimestamp* ts1=x1->GetTimestamp(); | |
1065 | if (!ts1) continue; | |
1066 | pos.Reset(); | |
1067 | time.Reset(); | |
1068 | x1->GetPosition(vec,"car"); | |
1069 | pos.Enter(vec[0],vec[1],vec[2]); | |
bfde1ea0 | 1070 | t0=fEvt->GetDifference(ts1,"ns"); |
25eefd00 | 1071 | time.Enter(t0); |
1072 | for (Int_t jet2=0; jet2<njets; jet2++) | |
67338e65 | 1073 | { |
25eefd00 | 1074 | jx2=(AliJet*)jets2.At(jet2); |
1075 | if (!jx2 || jet2==jet1) continue; | |
1076 | AliPosition* x2=jx2->GetReferencePoint(); | |
1077 | if (!x2) continue; | |
1078 | AliTimestamp* ts2=x2->GetTimestamp(); | |
1079 | if (!ts2) continue; | |
1080 | ang=jx1->GetOpeningAngle(*jx2,"deg"); | |
1081 | if (ang<=fJangmax) | |
67338e65 | 1082 | { |
25eefd00 | 1083 | if (!fJinvol) |
413d0114 | 1084 | { |
25eefd00 | 1085 | dist=jx1->GetDistance(jx2); |
1086 | } | |
1087 | else | |
1088 | { | |
1089 | dist=jx1->GetDistance(x2); | |
1090 | dist2=jx2->GetDistance(x1); | |
1091 | if (dist2<dist) dist=dist2; | |
1092 | } | |
1093 | if (dist<=fJdistmax) | |
1094 | { | |
1095 | x2->GetPosition(vec,"car"); | |
1096 | pos.Enter(vec[0],vec[1],vec[2]); | |
bfde1ea0 | 1097 | t0=fEvt->GetDifference(ts2,"ns"); |
25eefd00 | 1098 | time.Enter(t0); |
1099 | for (Int_t jtk=1; jtk<=jx2->GetNtracks(); jtk++) | |
1100 | { | |
bfde1ea0 | 1101 | AliTrack* te=jx2->GetTrack(jtk); |
25eefd00 | 1102 | if (!te) continue; |
1103 | jx1->AddTrack(te); | |
1104 | } | |
1105 | jets2.RemoveAt(jet2); | |
1106 | if (fJiterate) merged=1; | |
413d0114 | 1107 | } |
67338e65 | 1108 | } |
25eefd00 | 1109 | } // End of jet2 loop |
1110 | ||
1111 | // Set the reference point data for this jet | |
1112 | for (Int_t k=1; k<=3; k++) | |
1113 | { | |
1114 | vec[k-1]=pos.GetMean(k); | |
1115 | err[k-1]=pos.GetSigma(k); | |
67338e65 | 1116 | } |
25eefd00 | 1117 | r0.SetPosition(vec,"car"); |
1118 | r0.SetPositionErrors(err,"car"); | |
bfde1ea0 | 1119 | r0.SetTimestamp((AliTimestamp&)*fEvt); |
25eefd00 | 1120 | AliTimestamp* jt0=r0.GetTimestamp(); |
1121 | t0=time.GetMean(1); | |
1122 | jt0->Add(0,0,(int)t0); | |
1123 | jx1->SetReferencePoint(r0); | |
1124 | ||
1125 | nah=jx1->GetNsignals(); | |
1126 | if (nah>nahmax) nahmax=nah; | |
1127 | } // End of jet1 loop | |
1128 | ||
1129 | jets2.Compress(); | |
1130 | ||
bfde1ea0 | 1131 | // The sum of 0.15*(nah-nahmax) and average qtc value per track for each jet |
1132 | // will be stored as the jet energy to enable sorting on this value | |
25eefd00 | 1133 | for (Int_t jjet=0; jjet<njets; jjet++) |
1134 | { | |
1135 | AliJet* jx=(AliJet*)jets2.At(jjet); | |
1136 | if (!jx) continue; | |
1137 | nah=jx->GetNsignals(); | |
1138 | ntk=jx->GetNtracks(); | |
1139 | sortval=0.15*float(nah-nahmax); | |
1140 | if (ntk) sortval+=jx->GetMomentum()/float(ntk); | |
1141 | jx->SetScalar(sortval); | |
67338e65 | 1142 | } |
25eefd00 | 1143 | |
bfde1ea0 | 1144 | // Order the jets w.r.t. decreasing quality value |
1145 | TObjArray* ordered=fEvt->SortJets(-2,&jets2); | |
25eefd00 | 1146 | njets=ordered->GetEntries(); |
1147 | jets2.Clear(); | |
1148 | for (Int_t icopy=0; icopy<njets; icopy++) | |
413d0114 | 1149 | { |
25eefd00 | 1150 | jets2.Add(ordered->At(icopy)); |
413d0114 | 1151 | } |
25eefd00 | 1152 | } // End of iterative while loop |
67338e65 | 1153 | } |
bfde1ea0 | 1154 | } |
1155 | /////////////////////////////////////////////////////////////////////////// | |
1156 | void IceDwalk::StoreTracks(TObjArray& jets2) | |
1157 | { | |
67338e65 | 1158 | // Store every jet as a reconstructed track in the event structure. |
413d0114 | 1159 | // The jet 3-momentum (normalised to 1) and reference point |
1160 | // (i.e.the average r0 and t0 of the constituent tracks) will make up | |
1161 | // the final track parameters. | |
67338e65 | 1162 | // All the associated hits of all the constituent tracks of the jet |
1163 | // will be associated to the final track. | |
1164 | // In case the jet angular separation was set <0, only the jet with | |
1165 | // the maximum number of tracks (i.e. the first one in the array) | |
1166 | // will be used to form a track. This will allow comparison with | |
1167 | // the standard Sieglinde processing. | |
bfde1ea0 | 1168 | |
1169 | Int_t njets=jets2.GetEntries(); | |
1170 | ||
1171 | if (fTrackname=="") fTrackname=ClassName(); | |
1172 | TString title=ClassName(); | |
1173 | title+=" reco track"; | |
67338e65 | 1174 | AliTrack t; |
bfde1ea0 | 1175 | t.SetNameTitle(fTrackname.Data(),title.Data()); |
30672a96 | 1176 | t.SetCharge(fCharge); |
bfde1ea0 | 1177 | Ali3Vector p; |
67338e65 | 1178 | for (Int_t jet=0; jet<njets; jet++) |
1179 | { | |
1180 | AliJet* jx=(AliJet*)jets2.At(jet); | |
1181 | if (!jx) continue; | |
413d0114 | 1182 | AliPosition* ref=jx->GetReferencePoint(); |
1183 | if (!ref) continue; | |
bfde1ea0 | 1184 | fEvt->AddTrack(t); |
1185 | AliTrack* trk=fEvt->GetTrack(fEvt->GetNtracks()); | |
67338e65 | 1186 | if (!trk) continue; |
bfde1ea0 | 1187 | trk->SetId(fEvt->GetNtracks(1)+1); |
67338e65 | 1188 | p=jx->Get3Momentum(); |
1189 | p/=p.GetNorm(); | |
1190 | trk->Set3Momentum(p); | |
413d0114 | 1191 | trk->SetReferencePoint(*ref); |
1192 | AliTimestamp* tt0=ref->GetTimestamp(); | |
1193 | if (tt0) trk->SetTimestamp(*tt0); | |
67338e65 | 1194 | for (Int_t jt=1; jt<=jx->GetNtracks(); jt++) |
1195 | { | |
1196 | AliTrack* tx=jx->GetTrack(jt); | |
1197 | if (!tx) continue; | |
67338e65 | 1198 | for (Int_t is=1; is<=tx->GetNsignals(); is++) |
1199 | { | |
bfde1ea0 | 1200 | AliSignal* sx1=tx->GetSignal(is); |
d0120ca2 | 1201 | if (sx1) sx1->AddTrack(*trk); |
67338e65 | 1202 | } |
1203 | } | |
67338e65 | 1204 | |
25eefd00 | 1205 | // Only take the jet with the highest quality number |
67338e65 | 1206 | // (i.e. the first jet in the list) when the user had selected |
1207 | // this reconstruction mode. | |
413d0114 | 1208 | if (fJangmax<0) break; |
67338e65 | 1209 | } |
67338e65 | 1210 | } |
1211 | /////////////////////////////////////////////////////////////////////////// |