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