1 /*******************************************************************************
2 * Copyright(c) 2003, IceCube Experiment at the South Pole. All rights reserved.
4 * Author: The IceCube RALICE-based Offline Project.
5 * Contributors are mentioned in the code where appropriate.
7 * Permission to use, copy, modify and distribute this software and its
8 * documentation strictly for non-commercial purposes is hereby granted
9 * without fee, provided that the above copyright notice appears in all
10 * copies and that both the copyright notice and this permission notice
11 * appear in the supporting documentation.
12 * The authors make no claims about the suitability of this software for
13 * any purpose. It is provided "as is" without express or implied warranty.
14 *******************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////
20 // TTask derived class to perform (Amanda-like) direct walk track reconstruction.
21 // This class is kept to provide a procedure that can closely match the
22 // performance of the original Amanda direct walk-II algorithm as implemented
23 // in the Sieglinde reconstruction framework.
24 // For a more sophisticated direct walk reconstruction procedure, please refer
25 // to the class IceDwalk.
26 // The IceDwalkx procedure is based on the method described in the Amanda
27 // publication in Nuclear Instruments and Methods A524 (2004) 179-180.
28 // However, the Amanda method has been extended with the intention to
29 // take also multiple (muon) tracks within 1 event into account.
30 // This will not only provide a means to reconstruct muon bundles and
31 // multiple track events in IceCube, but will also allow to reduce the
32 // background of faked upgoing muons as a result of multiple downgoing
33 // muons hitting the top and bottom parts of the detector.
34 // A further extension of the original Amanda method is the separate treatment
35 // of the phase and group velocities as introduced in collaboration with
36 // George Japaridze (Clark Atlanta University, USA) which will provide more
37 // accurate time residuals due to the different velocities of the Cerenkov
38 // wave front (v_phase) and the actually detected photons (v_group).
39 // This distinction between v_phase and v_group can be (de)activated via the
40 // memberfunction SetVgroupUsage(). By default the distinction between v_phase
41 // and v_group is activated in the constructor of this class.
42 // To prevent waisting CPU time in trying to reconstruct (high-energy) cascade
43 // events, or to select specifically reconstruction of low multiplicity events,
44 // the user may invoke the memberfunctions SetMaxModA() and SetMinModA().
45 // This allows selection of events for processing with a certain maximum and/or
46 // minimum number of good Amanda OMs firing.
47 // By default the minimum and maximum are set to 0 and 999, respectively,
48 // in the constructor, which implies no multiplicity selection.
49 // The maximum number of good hits per Amanda OM to be used for the reconstruction
50 // can be specified via the memberfunction SetMaxHitsA().
51 // By default only the first good hit of each Amanda OM is used to be consistent
52 // with the original Sieglinde direct walk procedure of the above NIM article.
53 // Note that when all the good hits of an OM are used, this may lead to large
54 // processing time in case many noise and/or afterpulse signals are not
55 // recognised by the hit cleaning procedure.
57 // Information about the actual parameter settings can be found in the event
58 // structure itself via the device named "IceDwalkx".
60 // The various reconstruction steps are summarised as follows :
62 // 1) Construction of track elements (TE's).
63 // A track element is a straight line connecting two hits that
64 // appeared at some minimum distance d and within some maximum
65 // time difference dt.
66 // The default values for d and dt are given in eq. (20) of the
67 // NIM article, but can be modified by the appropriate Set functions.
68 // For dt a default margin of 30 ns is used (according to eq. (20)),
69 // but also this margin may be modified via the appropriate Set function.
70 // The reference point r0 of the TE is taken as the center between
71 // the two hit positions and the TE timestamp t0 at the position r0
72 // is taken as the IceEvent timestamp increased by the average of the
73 // two hit times. So, all timestamps contain the overall IceEvent
74 // timestamp as a basis. This means that time differences can be
75 // obtained via the AliTimestamp facilities (supporting upto picosecond
76 // precision when available).
77 // The TE direction is given by the relative position of the two hits.
79 // 2) Each TE will obtain so called associated hits.
80 // A hit is associated to a TE when it fulfills both the conditions
82 // -30 < tres < 300 ns
83 // dhit < 25*(tres+30)^(1/4) meter
85 // tres : time residual
86 // Difference between the observed hit time and the time expected
87 // for a direct photon hit.
88 // dhit : Distance between the hit and the TE
90 // 3) Construction of track candidates (TC's).
91 // These are TE's that fulfill both the conditions (see eq. (21) in the NIM article)
96 // qtc=min(nah,0.3*sigmal+7)
97 // qtcmax=max(qtc) of all TE's with sigmal>=20 meter
99 // where we have used the observables :
101 // nah : Number of associated hits.
102 // sigmal : rms variance of the distances between r0 and the projection
103 // points on the track of the various associated hit positions.
105 // Note : The following additional quality selection as indicated
106 // in the NIM article is not used anymore.
110 // 4) The remaining track candidates are clustered into jets when their directions
111 // are within a certain maximum opening angle.
112 // In addition the distance between their r0's must be below a certain maximum
113 // or the relative r0 direction must fall within a certain maximum opening angle
114 // w.r.t. the jet-starting track candidate.
115 // The latter criterion prevents clustering of (nearly) parallel track candidates
116 // crossing the detector a very different locations (e.g. muon bundles).
117 // The default maximum track opening angle is 15 degrees, but can be modified
118 // via the SetTangmax memberfunction.
119 // The remaining parameters related to the r0 criteria can be modified via
120 // the SetRtdmax and SetRtangmax memberfunctions.
121 // See the corresponding docs for defaults etc...
123 // The average of all the r0 and t0 values of the constituent TC's
124 // of the jet will provide the r0 and t0 (i.e. reference point) of the jet.
126 // 5) The jets are merged when their directions are within a certain maximum
127 // opening angle. Before the merging, all jets are ordered w.r.t. decreasing
128 // number of associated hits. This guarantees that the jet with the largest
129 // number of associated hits is the starting jet in the merging process.
130 // In addition the distance between their r0's must be below a certain maximum
131 // or the relative r0 direction must fall within a certain maximum opening angle
132 // w.r.t. the starting jet.
133 // The latter criterion prevents merging of (nearly) parallel tracks/jets
134 // crossing the detector a very different locations (e.g. muon bundles).
135 // The default maximum opening angle is half the TC maximum opening angle,
136 // but can be modified via the SetJangmax memberfunction.
137 // The remaining parameters related to the r0 criteria can be modified via
138 // the SetRjdmax and SetRjangmax memberfunctions.
139 // See the corresponding docs for defaults etc...
141 // Note : Setting the maximum jet opening angle to <=0 will prevent
142 // the merging of jets.
144 // The average of all the r0 and t0 values of the merged jets will provide
145 // the r0 and t0 (i.e. reference point) of the final jet.
147 // 6) The remaining jets are ordered w.r.t. decreasing number of associated hits.
148 // Each remaining jet will provide the parameters (e.g. direction)
149 // for a reconstructed track.
150 // The track 3-momentum is set to the total jet 3-momentum, normalised
151 // to 1 GeV. The mass and charge of the track are left 0.
152 // The reference point data of the jet will provide the r0 and t0
153 // (i.e. reference point) of the track.
155 // All these tracks will be stored in the IceEvent structure with as
156 // default "IceDwalkx" as the name of the track.
157 // This track name identifier can be modified by the user via the
158 // SetTrackName() memberfunction. This will allow unique identification
159 // of tracks which are produced when re-processing existing data with
160 // different criteria.
161 // By default the charge of the produced tracks is set to 0, since
162 // no distinction can be made between positive or negative tracks.
163 // However, the user can define the track charge by invokation
164 // of the memberfunction SetCharge().
165 // This facility may be used to distinguish tracks produced by the
166 // various reconstruction algorithms in a (3D) colour display
167 // (see the class AliHelix for further details).
169 // Note : In case the maximum jet opening angle was specified <0,
170 // only the jet with the maximum number of associated hits will
171 // appear as a reconstructed track in the IceEvent structure.
172 // This will allow comparison with the standard Sieglinde
173 // single track direct walk reconstruction results.
175 // For further details the user is referred to NIM A524 (2004) 169.
177 // Note : This algorithm works best on data which has been calibrated
178 // (IceCalibrate), cross talk corrected (IceXtalk) and cleaned
179 // from noise hits etc. (IceCleanHits).
181 //--- Author: Nick van Eijndhoven 07-oct-2005 Utrecht University
182 //- Modified: NvE $Date$ Utrecht University
183 ///////////////////////////////////////////////////////////////////////////
185 #include "IceDwalkx.h"
186 #include "Riostream.h"
188 ClassImp(IceDwalkx) // Class implementation to enable ROOT I/O
190 IceDwalkx::IceDwalkx(const char* name,const char* title) : TTask(name,title)
192 // Default constructor.
193 // The various reconstruction parameters are initialised to the values
194 // as mentioned in NIM A524 (2004) 179-180.
195 // The newly introduced angular separation parameter for jet merging
196 // is initialised as half the value of the angular separation parameter
197 // for track candidate clustering.
203 fJangmax=fTangmax/2.;
210 fTrackname="IceDwalkx";
213 ///////////////////////////////////////////////////////////////////////////
214 IceDwalkx::~IceDwalkx()
216 // Default destructor.
218 ///////////////////////////////////////////////////////////////////////////
219 void IceDwalkx::SetDmin(Float_t d)
221 // Set minimum hit distance (in m) to form a track element.
222 // In the constructor the default has been set to 50 meter, in accordance
223 // to eq.(20) of NIM A524 (2004) 179.
226 ///////////////////////////////////////////////////////////////////////////
227 void IceDwalkx::SetDtmarg(Int_t dt)
229 // Set maximum hit time difference margin (in ns) for track elements.
230 // In the constructor the default has been set to 30 ns, in accordance
231 // to eq.(20) of NIM A524 (2004) 179.
234 ///////////////////////////////////////////////////////////////////////////
235 void IceDwalkx::SetTangmax(Float_t ang)
237 // Set maximum angular separation (in deg) for track candidate clustering
239 // In the constructor the default has been set to 15 deg, in accordance
240 // to NIM A524 (2004) 180.
242 // Note : This function also sets automatically the value of the maximum
243 // angular separation for jet merging into 1 single track to ang/2.
244 // In order to specify a different max. jet merging separation angle,
245 // one has to invoke the memberfunction SetJangmax afterwards.
250 ///////////////////////////////////////////////////////////////////////////
251 void IceDwalkx::SetRtangmax(Float_t ang)
253 // Set maximum angular separation (in deg) for the relative direction of the
254 // two r0's of two track candidates (w.r.t. the direction of the starting
255 // track candidate) in the track clustering process.
256 // In the constructor the default has been set to 15 deg, corresponding
257 // to the default max. angular separation for track candidate clustering.
259 // Note : This function also sets automatically the value of the maximum
260 // angular separation for the relative direction of the two r0's
261 // of two jets (w.r.t. the direction of the starting jet)
262 // in the jet merging process.
263 // In order to specify a different value related to jet merging,
264 // one has to invoke the memberfunction SetRjangmax afterwards.
269 ///////////////////////////////////////////////////////////////////////////
270 void IceDwalkx::SetRtdmax(Float_t d)
272 // Set maximum distance (in m) of the two r0's of two track candidates
273 // in the track clustering process.
274 // This will allow clustering of tracks with very close r0's, of which
275 // their relative direction may point in any direction.
276 // In the constructor the default has been set 0 until further tuning
277 // of this parameter has been achieved.
279 // Note : In case the distance between the two r0's exceeds the maximum,
280 // the track candidates will still be clustered if the relative
281 // direction of the two r0's falls within the maximum separation
282 // angle w.r.t. the starting track direction.
286 ///////////////////////////////////////////////////////////////////////////
287 void IceDwalkx::SetJangmax(Float_t ang)
289 // Set angular separation (in deg) within which jets are merged into 1
291 // In the constructor the default has been set 7.5 deg, being half of the
292 // value of the default track candidate clustering separation angle.
296 // 1) Setting ang=0 will prevent jet merging.
297 // Consequently, every jet will appear as a separate track in the
298 // reconstruction result.
299 // 2) Setting ang<0 will prevent jet merging.
300 // In addition, only the jet with the maximum number of tracks will
301 // appear as a track in the reconstruction result.
302 // This situation resembles the standard Sieglinde direct walk processing
303 // and as such can be used to perform comparison studies.
307 ///////////////////////////////////////////////////////////////////////////
308 void IceDwalkx::SetRjangmax(Float_t ang)
310 // Set maximum angular separation (in deg) for the relative direction of the
311 // two r0's of two jets (w.r.t. the direction of the starting jet)
312 // in the jet merging process.
313 // In the constructor the default has been set to the corresponding value
314 // of the same parameter related to the track clustering.
318 ///////////////////////////////////////////////////////////////////////////
319 void IceDwalkx::SetRjdmax(Float_t d)
321 // Set maximum distance (in m) of the two r0's of two jets in the
322 // jet merging process.
323 // This will allow merging of jets with rather close r0's, of which
324 // their relative direction may point in any direction.
325 // In the constructor the default has been set to 50 meter, corresponding
326 // to the value of the minimum hit distance to form a track element.
328 // Note : In case the distance between the two r0's exceeds the maximum,
329 // the jets will still be merged if the relative direction of the
330 // two r0's falls within the maximum separation angle w.r.t. the
331 // starting jet direction.
335 ///////////////////////////////////////////////////////////////////////////
336 void IceDwalkx::SetMaxModA(Int_t nmax)
338 // Set the maximum number of good Amanda modules that may have fired
339 // in order to process this event.
340 // This allows suppression of processing (high-energy) cascade events
341 // with this direct walk tracking to prevent waisting cpu time for cases
342 // in which tracking doesn't make sense anyhow.
343 // Furthermore it allows selection of low multiplicity events for processing.
344 // By default the maximum number of Amanda modules is set to 999 in the ctor,
345 // which implies no selection on maximum module multiplicity.
346 // See also the memberfunction SetMinModA().
349 ///////////////////////////////////////////////////////////////////////////
350 void IceDwalkx::SetMinModA(Int_t nmin)
352 // Set the minimum number of good Amanda modules that must have fired
353 // in order to process this event.
354 // This allows selection of a minimal multiplicity for events to be processed.
355 // By default the minimum number of Amanda modules is set to 0 in the ctor,
356 // which implies no selection on minimum module multiplicity.
357 // See also the memberfunction SetMaxModA().
360 ///////////////////////////////////////////////////////////////////////////
361 void IceDwalkx::SetMaxHitsA(Int_t nmax)
363 // Set the maximum number of good hits per Amanda module to be processed.
366 // nmax = 0 : No maximum limit set; all available good hits will be processed
367 // < 0 : No hits will be processed
369 // In case the user selects a maximum number of good hits per module, all the
370 // hits of each module will be ordered w.r.t. increasing hit time (LE).
371 // This allows selection of processing e.g. only the first good hits etc...
372 // By default the maximum number of good hits per Amanda modules is set to 1
373 // in the ctor, which implies just processing only the first good hit of
377 ///////////////////////////////////////////////////////////////////////////
378 void IceDwalkx::SetVgroupUsage(Int_t flag)
380 // (De)activate the distinction between v_phase and v_group of the Cherenkov light.
382 // flag = 0 : No distinction between v_phase and v_group
383 // = 1 : Separate treatment of v_phase and v_group
385 // By default the distinction between v_phase and v_group is activated
386 // in the constructor of this class.
389 ///////////////////////////////////////////////////////////////////////////
390 void IceDwalkx::SetTrackName(TString s)
392 // Set (alternative) name identifier for the produced first guess tracks.
393 // This allows unique identification of (newly) produced direct walk tracks
394 // in case of re-processing of existing data with different criteria.
395 // By default the produced first guess tracks have the name "IceDwalkx"
396 // which is set in the constructor of this class.
399 ///////////////////////////////////////////////////////////////////////////
400 void IceDwalkx::SetCharge(Float_t charge)
402 // Set user defined charge for the produced first guess tracks.
403 // This allows identification of these tracks on color displays.
404 // By default the produced first guess tracks have charge=0
405 // which is set in the constructor of this class.
408 ///////////////////////////////////////////////////////////////////////////
409 void IceDwalkx::Exec(Option_t* opt)
411 // Implementation of the direct walk track reconstruction.
414 AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
418 IceEvent* evt=(IceEvent*)parent->GetObject("IceEvent");
421 // Enter the reco parameters as a device in the event
423 params.SetNameTitle("IceDwalkx","IceDwalkx reco parameters");
424 params.SetSlotName("Dmin",1);
425 params.SetSlotName("Dtmarg",2);
426 params.SetSlotName("Tangmax",3);
427 params.SetSlotName("Rtangmax",4);
428 params.SetSlotName("Rtdmax",5);
429 params.SetSlotName("Jangmax",6);
430 params.SetSlotName("Rjangmax",7);
431 params.SetSlotName("Rjdmax",8);
432 params.SetSlotName("MaxmodA",9);
433 params.SetSlotName("MinmodA",10);
434 params.SetSlotName("MaxhitsA",11);
435 params.SetSlotName("Vgroup",12);
437 params.SetSignal(fDmin,1);
438 params.SetSignal(fDtmarg,2);
439 params.SetSignal(fTangmax,3);
440 params.SetSignal(fRtangmax,4);
441 params.SetSignal(fRtdmax,5);
442 params.SetSignal(fJangmax,6);
443 params.SetSignal(fRjangmax,7);
444 params.SetSignal(fRjdmax,8);
445 params.SetSignal(fMaxmodA,9);
446 params.SetSignal(fMinmodA,10);
447 params.SetSignal(fMaxhitsA,11);
448 params.SetSignal(fVgroup,12);
450 evt->AddDevice(params);
452 if (fMaxhitsA<0) return;
454 // Fetch all fired Amanda OMs for this event
455 TObjArray* aoms=evt->GetDevices("IceAOM");
456 Int_t naoms=aoms->GetEntries();
459 // Check for the minimum and/or maximum number of good fired Amanda OMs
461 for (Int_t iom=0; iom<naoms; iom++)
463 IceGOM* omx=(IceGOM*)aoms->At(iom);
465 if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue;
468 if (ngood<fMinmodA || ngood>fMaxmodA) return;
470 const Float_t pi=acos(-1.);
471 const Float_t c=0.299792458; // Light speed in vacuum in meters per ns
472 const Float_t npice=1.31768387; // Phase refractive index (c/v_phase) of ice
473 const Float_t ngice=1.35075806; // Group refractive index (c/v_group) of ice
474 const Float_t thetac=acos(1./npice); // Cherenkov angle (in radians)
476 // Angular reduction of complement of thetac due to v_phase and v_group difference
478 if (fVgroup) alphac=atan((1.-npice/ngice)/sqrt(npice*npice-1.));
480 // Storage of track elements.
500 // Check the hits of Amanda OM pairs for posible track elements.
501 // Also all the good hits are stored in the meantime (to save CPU time)
502 // for hit association with the various track elements lateron.
504 for (Int_t i1=0; i1<naoms; i1++) // First OM of the pair
506 IceGOM* omx1=(IceGOM*)aoms->At(i1);
508 if (omx1->GetDeadValue("LE")) continue;
509 r1=omx1->GetPosition();
510 // Select all the good hits of this first OM
512 // Determine the max. number of hits to be processed for this OM
514 if (fMaxhitsA>0 && omx1->GetNhits()>fMaxhitsA) ordered=omx1->SortHits("LE",1,0,7);
516 for (Int_t j1=1; j1<=omx1->GetNhits(); j1++)
520 if (nh1>=fMaxhitsA) break;
521 sx1=(AliSignal*)ordered->At(j1-1);
525 sx1=omx1->GetHit(j1);
528 if (sx1->GetDeadValue("ADC") || sx1->GetDeadValue("LE") || sx1->GetDeadValue("TOT")) continue;
530 // Also store all good hits in the total hit array
535 // No further pair to be formed with the last OM in the list
536 if (i1==(naoms-1)) break;
538 nh1=hits1.GetEntries();
541 for (Int_t i2=i1+1; i2<naoms; i2++) // Second OM of the pair
543 IceGOM* omx2=(IceGOM*)aoms->At(i2);
545 if (omx2->GetDeadValue("LE")) continue;
546 r2=omx2->GetPosition();
550 if (dist<=fDmin) continue;
552 // Select all the good hits of this second OM
554 // Determine the max. number of hits to be processed for this OM
556 if (fMaxhitsA>0 && omx2->GetNhits()>fMaxhitsA) ordered=omx2->SortHits("LE",1,0,7);
558 for (Int_t j2=1; j2<=omx2->GetNhits(); j2++)
562 if (nh2>=fMaxhitsA) break;
563 sx2=(AliSignal*)ordered->At(j2-1);
567 sx2=omx2->GetHit(j2);
570 if (sx2->GetDeadValue("ADC") || sx2->GetDeadValue("LE") || sx2->GetDeadValue("TOT")) continue;
575 nh2=hits2.GetEntries();
578 // Position r0 in between the two OMs and normalised relative direction r12
580 r0.SetPosition((Ali3Vector&)rsum);
583 // Check all hit pair combinations of these two OMs for possible track elements
584 dtmax=dist/c+float(fDtmarg);
585 for (Int_t ih1=0; ih1<nh1; ih1++) // Hits of first OM
587 sx1=(AliSignal*)hits1.At(ih1);
589 for (Int_t ih2=0; ih2<nh2; ih2++) // Hits of second OM
591 sx2=(AliSignal*)hits2.At(ih2);
593 t1=sx1->GetSignal("LE",7);
594 t2=sx2->GetSignal("LE",7);
598 if (fabs(dt)>=dtmax) continue;
603 r0.SetTimestamp((AliTimestamp&)*evt);
604 AliTimestamp* tsx=r0.GetTimestamp();
605 tsx->Add(0,0,(int)t0);
606 te->SetReferencePoint(r0);
607 te->Set3Momentum(r12);
610 } // end of loop over the second OM of the pair
611 } // end of loop over first OM of the pair
613 // Association of hits to the various track elements
614 // For the time being all track elements will be treated,
615 // but in a later stage one could select only the TE's of a certain
616 // 3 ns margin slot in the TE map to save CPU time.
617 Int_t nte=tes.GetEntries();
618 Int_t nh=hits.GetEntries();
622 AliSample levers; // Statistics of the assoc. hit lever arms
623 AliSignal fit; // Storage of Q value etc... for each track candidate
624 fit.SetSlotName("QTC",1);
625 fit.SetSlotName("SIGMAL",2);
626 Float_t qtc=0,qmax=0;
627 Int_t nah; // Number of associated hits for a certain TE
628 Float_t sigmal; // The mean lever arm of the various associated hits
629 Float_t hproj; // Projected hit position on the track w.r.t. r0
630 for (Int_t jte=0; jte<nte; jte++)
632 te=(AliTrack*)tes.At(jte);
634 AliPosition* tr0=te->GetReferencePoint();
635 AliTimestamp* tt0=tr0->GetTimestamp();
636 t0=evt->GetDifference(tt0,"ns");
637 p=te->Get3Momentum();
639 for (Int_t jh=0; jh<nh; jh++)
641 sx1=(AliSignal*)hits.At(jh);
643 IceGOM* omx=(IceGOM*)sx1->GetDevice();
645 r1=omx->GetPosition();
646 d=te->GetDistance(r1);
649 dist=hproj+d/tan(pi/2.-thetac-alphac);
651 t1=sx1->GetSignal("LE",7);
654 if (tres<-30 || tres>300 || d>25.*pow(tres+30.,0.25)) continue;
656 // Associate this hit to the TE
658 levers.Enter(fabs(hproj));
661 // Determine the Q quality of the various TE's.
662 // Good quality TE's will be called track candidates (TC's)
663 nah=te->GetNsignals();
664 sigmal=levers.GetSigma(1);
666 if (qtc>nah) qtc=nah;
667 if (sigmal<20) qtc=-1; // Reject TE's with sigmal<20 meter
668 fit.SetSignal(qtc,"QTC");
669 fit.SetSignal(sigmal,"SIGMAL");
670 te->SetFitDetails(fit);
671 if (qtc>qmax) qmax=qtc;
674 // Perform selection on Q value in case of multiple track candidates
675 for (Int_t jtc=0; jtc<nte; jtc++)
677 te=(AliTrack*)tes.At(jtc);
679 nah=te->GetNsignals();
680 sx1=(AliSignal*)te->GetFitDetails();
685 qtc=sx1->GetSignal("QTC");
686 sigmal=sx1->GetSignal("SIGMAL");
695 nte=tes.GetEntries();
699 // Cluster track candidates within a certain opening angle into jets.
700 // Also the relative direction of the both r0's of the track candidates
701 // should be within a certain opening angle w.r.t. the starting track direction,
702 // unless the distance between the two r0's is below a certain maximum.
703 // The latter prevents clustering of (nearly) parallel track candidates
704 // crossing the detector a very different locations (e.g. muon bundles).
705 // The average r0 and t0 of the constituent tracks will be taken as the
706 // jet reference point.
713 Float_t vec[3],err[3];
714 for (Int_t jtc1=0; jtc1<nte; jtc1++)
716 te=(AliTrack*)tes.At(jtc1);
718 AliPosition* x1=te->GetReferencePoint();
720 AliTimestamp* ts1=x1->GetTimestamp();
722 AliJet* jx=new AliJet();
726 x1->GetPosition(vec,"car");
727 pos.Enter(vec[0],vec[1],vec[2]);
728 t0=evt->GetDifference(ts1,"ns");
730 for (Int_t jtc2=0; jtc2<nte; jtc2++)
732 if (jtc2==jtc1) continue;
733 te2=(AliTrack*)tes.At(jtc2);
735 ang=te->GetOpeningAngle(*te2,"deg");
738 AliPosition* x2=te2->GetReferencePoint();
740 AliTimestamp* ts2=x2->GetTimestamp();
742 dist=x1->GetDistance(x2);
743 dt=ts1->GetDifference(ts2,"ns");
748 ang=te->GetOpeningAngle(r12,"deg");
749 if (ang<=fRtangmax || dist<fRtdmax)
751 x2->GetPosition(vec,"car");
752 pos.Enter(vec[0],vec[1],vec[2]);
753 t0=evt->GetDifference(ts2,"ns");
761 // Set the reference point data for this jet
762 for (Int_t j=1; j<=3; j++)
764 vec[j-1]=pos.GetMean(j);
765 err[j-1]=pos.GetSigma(j);
767 r0.SetPosition(vec,"car");
768 r0.SetPositionErrors(err,"car");
769 r0.SetTimestamp((AliTimestamp&)*evt);
770 AliTimestamp* jt0=r0.GetTimestamp();
772 jt0->Add(0,0,(int)t0);
773 jx->SetReferencePoint(r0);
775 // Store this jet for further processing if ntracks>1
776 if (jx->GetNtracks() > 1 || fTangmax<=0 || fRtangmax<=0)
780 else // Only keep single-track jets which have qtc=qmax
782 sx1=(AliSignal*)te->GetFitDetails();
784 if (sx1) qtc=sx1->GetSignal("QTC");
785 if (qtc>=(qmax-1.e-10))
795 Int_t njets=jets.GetEntries();
799 // Order the jets w.r.t. decreasing number of associated hits
800 ordered=evt->SortJets(-12,&jets);
801 TObjArray jets2(*ordered);
803 // Merge jets within a certain opening to provide the final track(s).
809 for (Int_t jet1=0; jet1<njets; jet1++)
811 jx1=(AliJet*)jets2.At(jet1);
813 AliPosition* x1=jx1->GetReferencePoint();
815 AliTimestamp* ts1=x1->GetTimestamp();
819 x1->GetPosition(vec,"car");
820 pos.Enter(vec[0],vec[1],vec[2]);
821 t0=evt->GetDifference(ts1,"ns");
823 for (Int_t jet2=jet1+1; jet2<njets; jet2++)
825 jx2=(AliJet*)jets2.At(jet2);
827 AliPosition* x2=jx2->GetReferencePoint();
829 AliTimestamp* ts2=x2->GetTimestamp();
831 ang=jx1->GetOpeningAngle(*jx2,"deg");
834 dist=x1->GetDistance(x2);
835 edist=x1->GetResultError();
836 dt=ts1->GetDifference(ts2,"ns");
839 ang=jx1->GetOpeningAngle(r12,"deg");
840 if (ang<=fRjangmax || dist<=(fRjdmax+edist))
842 x2->GetPosition(vec,"car");
843 pos.Enter(vec[0],vec[1],vec[2]);
844 t0=evt->GetDifference(ts2,"ns");
846 for (Int_t jtk=1; jtk<=jx2->GetNtracks(); jtk++)
848 te=jx2->GetTrack(jtk);
852 jets2.RemoveAt(jet2);
856 // Set the reference point data for this jet
857 for (Int_t k=1; k<=3; k++)
859 vec[k-1]=pos.GetMean(k);
860 err[k-1]=pos.GetSigma(k);
862 r0.SetPosition(vec,"car");
863 r0.SetPositionErrors(err,"car");
864 r0.SetTimestamp((AliTimestamp&)*evt);
865 AliTimestamp* jt0=r0.GetTimestamp();
867 jt0->Add(0,0,(int)t0);
868 jx1->SetReferencePoint(r0);
871 njets=jets2.GetEntries();
874 // Order the merged jets w.r.t. decreasing number of associated hits
875 ordered=evt->SortJets(-12,&jets2);
876 TObjArray jets3(*ordered);
878 // Store every jet as a reconstructed track in the event structure.
879 // The jet 3-momentum (normalised to 1) and reference point
880 // (i.e.the average r0 and t0 of the constituent tracks) will make up
881 // the final track parameters.
882 // All the associated hits of all the constituent tracks of the jet
883 // will be associated to the final track.
884 // In case the jet angular separation was set <0, only the jet with
885 // the maximum number of tracks (i.e. the first one in the array)
886 // will be used to form a track. This will allow comparison with
887 // the standard Sieglinde processing.
889 t.SetNameTitle(fTrackname.Data(),"IceDwalkx direct walk track");
890 t.SetCharge(fCharge);
891 for (Int_t jet=0; jet<njets; jet++)
893 AliJet* jx=(AliJet*)jets3.At(jet);
895 AliPosition* ref=jx->GetReferencePoint();
898 AliTrack* trk=evt->GetTrack(evt->GetNtracks());
900 trk->SetId(evt->GetNtracks(1)+1);
901 p=jx->Get3Momentum();
903 trk->Set3Momentum(p);
904 trk->SetReferencePoint(*ref);
905 AliTimestamp* tt0=ref->GetTimestamp();
906 if (tt0) trk->SetTimestamp(*tt0);
907 for (Int_t jt=1; jt<=jx->GetNtracks(); jt++)
909 AliTrack* tx=jx->GetTrack(jt);
911 for (Int_t is=1; is<=tx->GetNsignals(); is++)
913 sx1=tx->GetSignal(is);
914 if (sx1) sx1->AddTrack(*trk);
918 // Only take the jet with the maximum number of associated hits
919 // (i.e. the first jet in the list) when the user had selected
920 // this reconstruction mode.
921 if (fJangmax<0) break;
924 ///////////////////////////////////////////////////////////////////////////