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 direct walk track reconstruction.
22 // In case an event has been rejected by an AliEventSelector (based) processor,
23 // this task (and its sub-tasks) is not executed.
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
32 // muons hitting the top and bottom parts of the detector.
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.
41 // To prevent waisting CPU time in trying to reconstruct (high-energy) cascade
42 // events, or to select specifically reconstruction of low multiplicity events,
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,
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
49 // can be specified via the memberfunction SetMaxHitsA().
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
52 // processing time in case many noise and/or afterpulse signals are not
53 // recognised by the hit cleaning procedure.
55 // Information about the actual parameter settings can be found in the event
56 // structure itself via the device named "IceDwalk".
58 // The various reconstruction steps are summarised as follows :
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
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().
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.
77 // 2) Each TE will obtain so called associated hits.
78 // A hit is associated to a TE when it fulfills both the conditions
80 // -30 < tres < 300 ns
83 // tres : time residual
84 // Difference between the observed hit time and the time expected
85 // for a direct photon hit.
86 // dhit : Distance traveled by the cherenkov photon from the track to the hit position
87 // lambda : Photon scattering length in ice
89 // By default F is set to 3.07126, but this can be modified via the memberfunction
92 // 3) Construction of track candidates (TC's).
93 // These are TE's that fulfill both the conditions
98 // where we have defined :
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.
104 // The track quality number qtc is defined as follows :
106 // qtc=nah*(term1+term2)-term3-term4-term5
108 // here we have defined :
110 // term1=2*spread/span
111 // term2=2*spreadL/spanL
112 // term3=|spread-expspread|/spread
113 // term4=|spreadL-expspreadL|/spreadL
114 // term5=|medianT|/spreadT
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|.
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| >
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.
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).
140 // 4) The remaining track candidates are clustered into jets when their directions
141 // are within a certain maximum opening angle.
142 // In addition a track candidate must within a certain maximum distance
143 // of the jet starting TC in order to get clustered.
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.
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.
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.
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 :
163 // qtcjet=|jet momentum|/ntracks
165 // This jet quality number is refined on basis of the number of hits
166 // associated to the jet as :
168 // qtcjet=qtcjet+0.2*(nah-nahmax)
170 // where we have defined :
172 // nah : Number of associated hits for the specific jet.
173 // nahmax : Maximum number of associated hits encountered for the jets.
175 // This qtcjet value is then used to order the various jets w.r.t.
176 // decreasing qtcjet quality number.
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.
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.
187 // The latter criterion prevents merging of (nearly) parallel tracks/jets
188 // crossing the detector a very different locations (e.g. muon bundles).
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.
193 // The default maximum opening angle is half the TC maximum opening angle,
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.
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.
204 // Note : Setting the maximum jet opening angle to <=0 will prevent
205 // the merging of jets.
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.
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.
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.
218 // The reference point data of the jet will provide the r0 and t0
219 // (i.e. reference point) of the track.
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.
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).
235 // Note : In case the maximum jet opening angle was specified <0,
236 // only the jet with the highest quality number will appear
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.
241 // For further details the user is referred to NIM A524 (2004) 169.
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).
247 //--- Author: Nick van Eijndhoven 07-oct-2005 Utrecht University
248 //- Modified: NvE $Date$ Utrecht University
249 ///////////////////////////////////////////////////////////////////////////
251 #include "IceDwalk.h"
252 #include "Riostream.h"
254 ClassImp(IceDwalk) // Class implementation to enable ROOT I/O
256 IceDwalk::IceDwalk(const char* name,const char* title) : TTask(name,title)
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.
271 fJangmax=fTangmax/2.;
282 ///////////////////////////////////////////////////////////////////////////
283 IceDwalk::~IceDwalk()
285 // Default destructor.
287 ///////////////////////////////////////////////////////////////////////////
288 void IceDwalk::SetDmin(Float_t d)
290 // Set minimum hit distance (in m) to form a track element.
291 // In the constructor the default has been set to 75 meter.
294 ///////////////////////////////////////////////////////////////////////////
295 void IceDwalk::SetDtmarg(Int_t dt)
297 // Set maximum hit time difference margin (in ns) for track elements.
298 // In the constructor the default has been set to 0 ns.
301 ///////////////////////////////////////////////////////////////////////////
302 void IceDwalk::SetMaxDhit(Float_t d)
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.
308 ///////////////////////////////////////////////////////////////////////////
309 void IceDwalk::SetTangmax(Float_t ang)
311 // Set maximum angular separation (in deg) for track candidate clustering
313 // In the constructor the default has been set to 15 deg, in accordance
314 // to NIM A524 (2004) 180.
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.
324 ///////////////////////////////////////////////////////////////////////////
325 void IceDwalk::SetTdistmax(Float_t d,Int_t invol)
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.
342 ///////////////////////////////////////////////////////////////////////////
343 void IceDwalk::SetJangmax(Float_t ang,Int_t iter)
345 // Set angular separation (in deg) within which jets are merged into 1
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.
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.
371 ///////////////////////////////////////////////////////////////////////////
372 void IceDwalk::SetJdistmax(Float_t d,Int_t invol)
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.
388 ///////////////////////////////////////////////////////////////////////////
389 void IceDwalk::SetMaxModA(Int_t nmax)
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,
398 // which implies no selection on maximum module multiplicity.
399 // See also the memberfunction SetMinModA().
402 ///////////////////////////////////////////////////////////////////////////
403 void IceDwalk::SetMinModA(Int_t nmin)
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().
413 ///////////////////////////////////////////////////////////////////////////
414 void IceDwalk::SetMaxHitsA(Int_t nmax)
416 // Set the maximum number of good hits per Amanda module to be processed.
419 // nmax = 0 : No maximum limit set; all good hits will be processed
420 // < 0 : No hits will be processed
422 // In case the user selects a maximum number of good hits per module, all the
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...
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.
429 ///////////////////////////////////////////////////////////////////////////
430 void IceDwalk::SetVgroupUsage(Int_t flag)
432 // (De)activate the distinction between v_phase and v_group of the Cherenkov light.
434 // flag = 0 : No distinction between v_phase and v_group
435 // = 1 : Separate treatment of v_phase and v_group
437 // By default the distinction between v_phase and v_group is activated
438 // in the constructor of this class.
441 ///////////////////////////////////////////////////////////////////////////
442 void IceDwalk::SetTrackName(TString s)
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.
447 // By default the produced first guess tracks have the name of the class
448 // by which they were produced.
451 ///////////////////////////////////////////////////////////////////////////
452 void IceDwalk::SetCharge(Float_t charge)
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.
460 ///////////////////////////////////////////////////////////////////////////
461 void IceDwalk::Exec(Option_t* opt)
463 // Implementation of the direct walk track reconstruction.
466 AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
470 fEvt=(IceEvent*)parent->GetObject("IceEvent");
473 // Only process accepted events
474 AliDevice* seldev=(AliDevice*)fEvt->GetDevice("AliEventSelector");
477 if (seldev->GetSignal("Select") < 0.1) return;
480 // Enter the reco parameters as a device in the event
482 params.SetNameTitle(ClassName(),"Reco parameters");
483 params.SetSlotName("Dmin",1);
484 params.SetSlotName("Dtmarg",2);
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);
498 params.SetSignal(fDmin,1);
499 params.SetSignal(fDtmarg,2);
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);
513 fEvt->AddDevice(params);
515 if (fMaxhitsA<0) return;
517 // Fetch all fired Amanda OMs for this event
518 TObjArray* aoms=fEvt->GetDevices("IceAOM");
520 Int_t naoms=aoms->GetEntries();
523 // Check for the minimum and/or maximum number of good fired Amanda OMs
525 for (Int_t iom=0; iom<naoms; iom++)
527 IceGOM* omx=(IceGOM*)aoms->At(iom);
529 if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue;
532 if (ngood<fMinmodA || ngood>fMaxmodA) return;
534 const Float_t c=0.299792458; // Light speed in vacuum in meters per ns
536 // Storage of track elements.
556 // Check the hits of Amanda OM pairs for possible track elements.
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.
560 for (Int_t i1=0; i1<naoms; i1++) // First OM of the pair
562 IceGOM* omx1=(IceGOM*)aoms->At(i1);
564 if (omx1->GetDeadValue("LE")) continue;
565 r1=omx1->GetPosition();
566 // Select all the good hits of this first OM
568 // Determine the max. number of hits to be processed for this OM
570 if (fMaxhitsA>0 && omx1->GetNhits()>fMaxhitsA) ordered=omx1->SortHits("LE",1,0,7);
572 for (Int_t j1=1; j1<=omx1->GetNhits(); j1++)
576 if (nh1>=fMaxhitsA) break;
577 sx1=(AliSignal*)ordered->At(j1-1);
581 sx1=omx1->GetHit(j1);
584 if (sx1->GetDeadValue("ADC") || sx1->GetDeadValue("LE") || sx1->GetDeadValue("TOT")) continue;
586 // Also store all good hits in the total hit array
591 // No further pair to be formed with the last OM in the list
592 if (i1==(naoms-1)) break;
594 nh1=hits1.GetEntries();
597 for (Int_t i2=i1+1; i2<naoms; i2++) // Second OM of the pair
599 IceGOM* omx2=(IceGOM*)aoms->At(i2);
601 if (omx2->GetDeadValue("LE")) continue;
602 r2=omx2->GetPosition();
606 if (dist<=fDmin) continue;
608 // Select all the good hits of this second OM
610 // Determine the max. number of hits to be processed for this OM
612 if (fMaxhitsA>0 && omx2->GetNhits()>fMaxhitsA) ordered=omx2->SortHits("LE",1,0,7);
614 for (Int_t j2=1; j2<=omx2->GetNhits(); j2++)
618 if (nh2>=fMaxhitsA) break;
619 sx2=(AliSignal*)ordered->At(j2-1);
623 sx2=omx2->GetHit(j2);
626 if (sx2->GetDeadValue("ADC") || sx2->GetDeadValue("LE") || sx2->GetDeadValue("TOT")) continue;
631 nh2=hits2.GetEntries();
634 // Position r0 in between the two OMs and normalised relative direction r12
636 r0.SetPosition((Ali3Vector&)rsum);
639 // Check all hit pair combinations of these two OMs for possible track elements
640 dtmax=dist/c+float(fDtmarg);
641 for (Int_t ih1=0; ih1<nh1; ih1++) // Hits of first OM
643 sx1=(AliSignal*)hits1.At(ih1);
645 for (Int_t ih2=0; ih2<nh2; ih2++) // Hits of second OM
647 sx2=(AliSignal*)hits2.At(ih2);
649 t1=sx1->GetSignal("LE",7);
650 t2=sx2->GetSignal("LE",7);
654 if (fabs(dt)>=dtmax) continue;
659 r0.SetTimestamp((AliTimestamp&)*fEvt);
660 AliTimestamp* tsx=r0.GetTimestamp();
661 tsx->Add(0,0,(int)t0);
662 te->SetReferencePoint(r0);
663 te->Set3Momentum(r12);
666 } // end of loop over the second OM of the pair
667 } // end of loop over first OM of the pair
669 // Association of hits to the various track elements.
672 AssociateHits(tes,hits,qmax,nahmax);
674 // Selection on quality (Q value) in case of multiple track candidates
675 SelectQvalue(tes,qmax);
677 Int_t nte=tes.GetEntries();
680 // Clustering of track candidates into jets
683 ClusterTracks(tes,jets,qmax);
685 Int_t njets=jets.GetEntries();
688 // Order the jets w.r.t. decreasing quality value
689 ordered=fEvt->SortJets(-2,&jets);
690 TObjArray jets2(*ordered);
695 // Production and storage of the final tracks
698 ///////////////////////////////////////////////////////////////////////////
699 void IceDwalk::AssociateHits(TObjArray& tes,TObjArray& hits,Float_t& qmax,Int_t& nahmax)
701 // Association of hits to the various track elements.
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)
710 // Angular reduction of complement of thetac due to v_phase and v_group difference
712 if (fVgroup) alphac=atan((1.-npice/ngice)/sqrt(npice*npice-1.));
714 Int_t nte=tes.GetEntries();
715 Int_t nh=hits.GetEntries();
722 Float_t dist,t0,tgeo,tres;
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");
753 Int_t nah; // Number of associated hits for a certain TE
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;
760 for (Int_t jte=0; jte<nte; jte++)
762 AliTrack* te=(AliTrack*)tes.At(jte);
764 AliPosition* tr0=te->GetReferencePoint();
765 AliTimestamp* tt0=tr0->GetTimestamp();
766 t0=fEvt->GetDifference(tt0,"ns");
767 p=te->Get3Momentum();
771 for (Int_t jh=0; jh<nh; jh++)
773 AliSignal* sx1=(AliSignal*)hits.At(jh);
775 IceGOM* omx=(IceGOM*)sx1->GetDevice();
777 r1=omx->GetPosition();
778 d=te->GetDistance(r1);
781 dist=hproj+d/tan(pi/2.-thetac-alphac);
783 t1=sx1->GetSignal("LE",7);
786 d=d/sin(thetac); // The distance traveled by a cherenkov photon
788 if (tres<-30 || tres>300 || d>fMaxdhit*lambda) continue;
790 // Associate this hit to the TE
792 levers.Enter(fabs(hproj));
797 // Determine the Q quality of the various TE's.
798 // Good quality TE's will be called track candidates (TC's)
799 nah=te->GetNsignals();
800 if (nah>nahmax) nahmax=nah;
801 lmin=levers.GetMinimum(1);
802 lmax=levers.GetMaximum(1);
804 medianl=levers.GetMedian(1);
805 meanl=levers.GetMean(1);
806 sigmal=levers.GetSigma(1);
807 spreadl=levers.GetSpread(1);
808 // Expected spread for a flat distribution
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
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);
827 if (span>0) term1=2.*spread/span;
830 if (spanl>0) term2=2.*spreadl/spanl;
833 if (spread>0) term3=fabs(spread-expspread)/spread;
836 if (spreadl>0) term4=fabs(spreadl-expspreadl)/spreadl;
839 if (spreadt>0) term5=fabs(mediant)/spreadt;
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
844 fit.SetSignal(qtc,"QTC");
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");
866 te->SetFitDetails(fit);
867 if (qtc>qmax) qmax=qtc;
870 ///////////////////////////////////////////////////////////////////////////
871 void IceDwalk::SelectQvalue(TObjArray& tes,Float_t qmax)
873 // Perform selection on Q value in case of multiple track candidates
875 Int_t nte=tes.GetEntries();
879 for (Int_t jtc=0; jtc<nte; jtc++)
881 AliTrack* te=(AliTrack*)tes.At(jtc);
883 nah=te->GetNsignals();
884 AliSignal* sx1=(AliSignal*)te->GetFitDetails();
886 if (sx1) qtc=sx1->GetSignal("QTC");
887 if (!nah || qtc<0.8*qmax)
892 else // Set Q value as momentum to provide a weight for jet clustering
896 p=te->Get3Momentum();
904 ///////////////////////////////////////////////////////////////////////////
905 void IceDwalk::ClusterTracks(TObjArray& tes,TObjArray& jets,Float_t qmax)
907 // Cluster track candidates within a certain opening angle into jets.
908 // Also the track should be within a certain maximum distance of the
909 // starting track in order to get clustered.
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.
915 Int_t nte=tes.GetEntries();
919 Float_t vec[3],err[3];
921 Float_t t0,dist,dist2;
922 Int_t nah=0,nahmax=0; // Determine the max. number of associated hits for the jets
924 for (Int_t jtc1=0; jtc1<nte; jtc1++)
926 AliTrack* te=(AliTrack*)tes.At(jtc1);
928 AliPosition* x1=te->GetReferencePoint();
930 AliTimestamp* ts1=x1->GetTimestamp();
932 AliJet* jx=new AliJet();
936 x1->GetPosition(vec,"car");
937 pos.Enter(vec[0],vec[1],vec[2]);
938 t0=fEvt->GetDifference(ts1,"ns");
940 for (Int_t jtc2=0; jtc2<nte; jtc2++)
942 if (jtc2==jtc1) continue;
943 AliTrack* te2=(AliTrack*)tes.At(jtc2);
945 ang=te->GetOpeningAngle(*te2,"deg");
948 AliPosition* x2=te2->GetReferencePoint();
950 AliTimestamp* ts2=x2->GetTimestamp();
954 dist=te->GetDistance(te2);
958 dist=te->GetDistance(x2);
959 dist2=te2->GetDistance(x1);
960 if (dist2<dist) dist=dist2;
964 x2->GetPosition(vec,"car");
965 pos.Enter(vec[0],vec[1],vec[2]);
966 t0=fEvt->GetDifference(ts2,"ns");
973 // Set the reference point data for this jet
974 for (Int_t j=1; j<=3; j++)
976 vec[j-1]=pos.GetMean(j);
977 err[j-1]=pos.GetSigma(j);
979 r0.SetPosition(vec,"car");
980 r0.SetPositionErrors(err,"car");
981 r0.SetTimestamp((AliTimestamp&)*fEvt);
982 AliTimestamp* jt0=r0.GetTimestamp();
984 jt0->Add(0,0,(int)t0);
985 jx->SetReferencePoint(r0);
987 // Store this jet for further processing if ntracks>1
988 if (jx->GetNtracks() > 1 || fTangmax<=0)
991 nah=jx->GetNsignals();
992 if (nah>nahmax) nahmax=nah;
994 else // Only keep single-track jets which have qtc=qmax
996 AliSignal* sx1=(AliSignal*)te->GetFitDetails();
998 if (sx1) qtc=sx1->GetSignal("QTC");
999 if (qtc>=(qmax-1.e-10))
1002 nah=jx->GetNsignals();
1003 if (nah>nahmax) nahmax=nah;
1012 Int_t njets=jets.GetEntries();
1015 // The sum of 0.15*(nah-nahmax) and average qtc value per track for each jet
1016 // will be stored as the jet energy to enable sorting on this value lateron
1019 for (Int_t ijet=0; ijet<njets; ijet++)
1021 AliJet* jx=(AliJet*)jets.At(ijet);
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);
1030 ///////////////////////////////////////////////////////////////////////////
1031 void IceDwalk::MergeJets(TObjArray& jets2)
1033 // Merge jets within a certain opening to provide the final track(s).
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.
1041 Int_t njets=jets2.GetEntries();
1045 Int_t ntk,nah,nahmax;
1046 Float_t ang,dist,dist2,t0;
1050 Float_t vec[3],err[3];
1058 for (Int_t jet1=0; jet1<njets; jet1++)
1060 jx1=(AliJet*)jets2.At(jet1);
1062 AliPosition* x1=jx1->GetReferencePoint();
1064 AliTimestamp* ts1=x1->GetTimestamp();
1068 x1->GetPosition(vec,"car");
1069 pos.Enter(vec[0],vec[1],vec[2]);
1070 t0=fEvt->GetDifference(ts1,"ns");
1072 for (Int_t jet2=0; jet2<njets; jet2++)
1074 jx2=(AliJet*)jets2.At(jet2);
1075 if (!jx2 || jet2==jet1) continue;
1076 AliPosition* x2=jx2->GetReferencePoint();
1078 AliTimestamp* ts2=x2->GetTimestamp();
1080 ang=jx1->GetOpeningAngle(*jx2,"deg");
1085 dist=jx1->GetDistance(jx2);
1089 dist=jx1->GetDistance(x2);
1090 dist2=jx2->GetDistance(x1);
1091 if (dist2<dist) dist=dist2;
1093 if (dist<=fJdistmax)
1095 x2->GetPosition(vec,"car");
1096 pos.Enter(vec[0],vec[1],vec[2]);
1097 t0=fEvt->GetDifference(ts2,"ns");
1099 for (Int_t jtk=1; jtk<=jx2->GetNtracks(); jtk++)
1101 AliTrack* te=jx2->GetTrack(jtk);
1105 jets2.RemoveAt(jet2);
1106 if (fJiterate) merged=1;
1109 } // End of jet2 loop
1111 // Set the reference point data for this jet
1112 for (Int_t k=1; k<=3; k++)
1114 vec[k-1]=pos.GetMean(k);
1115 err[k-1]=pos.GetSigma(k);
1117 r0.SetPosition(vec,"car");
1118 r0.SetPositionErrors(err,"car");
1119 r0.SetTimestamp((AliTimestamp&)*fEvt);
1120 AliTimestamp* jt0=r0.GetTimestamp();
1122 jt0->Add(0,0,(int)t0);
1123 jx1->SetReferencePoint(r0);
1125 nah=jx1->GetNsignals();
1126 if (nah>nahmax) nahmax=nah;
1127 } // End of jet1 loop
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
1133 for (Int_t jjet=0; jjet<njets; jjet++)
1135 AliJet* jx=(AliJet*)jets2.At(jjet);
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);
1144 // Order the jets w.r.t. decreasing quality value
1145 TObjArray* ordered=fEvt->SortJets(-2,&jets2);
1146 njets=ordered->GetEntries();
1148 for (Int_t icopy=0; icopy<njets; icopy++)
1150 jets2.Add(ordered->At(icopy));
1152 } // End of iterative while loop
1155 ///////////////////////////////////////////////////////////////////////////
1156 void IceDwalk::StoreTracks(TObjArray& jets2)
1158 // Store every jet as a reconstructed track in the event structure.
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.
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.
1169 Int_t njets=jets2.GetEntries();
1171 if (fTrackname=="") fTrackname=ClassName();
1172 TString title=ClassName();
1173 title+=" reco track";
1175 t.SetNameTitle(fTrackname.Data(),title.Data());
1176 t.SetCharge(fCharge);
1178 for (Int_t jet=0; jet<njets; jet++)
1180 AliJet* jx=(AliJet*)jets2.At(jet);
1182 AliPosition* ref=jx->GetReferencePoint();
1185 AliTrack* trk=fEvt->GetTrack(fEvt->GetNtracks());
1187 trk->SetId(fEvt->GetNtracks(1)+1);
1188 p=jx->Get3Momentum();
1190 trk->Set3Momentum(p);
1191 trk->SetReferencePoint(*ref);
1192 AliTimestamp* tt0=ref->GetTimestamp();
1193 if (tt0) trk->SetTimestamp(*tt0);
1194 for (Int_t jt=1; jt<=jx->GetNtracks(); jt++)
1196 AliTrack* tx=jx->GetTrack(jt);
1198 for (Int_t is=1; is<=tx->GetNsignals(); is++)
1200 AliSignal* sx1=tx->GetSignal(is);
1201 if (sx1) sx1->AddTrack(*trk);
1205 // Only take the jet with the highest quality number
1206 // (i.e. the first jet in the list) when the user had selected
1207 // this reconstruction mode.
1208 if (fJangmax<0) break;
1211 ///////////////////////////////////////////////////////////////////////////