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.
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
28 // muons hitting the top and bottom parts of the detector.
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.
37 // To prevent waisting CPU time in trying to reconstruct (high-energy) cascade
38 // events, or to select specifically reconstruction of low multiplicity events,
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,
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
45 // can be specified via the memberfunction SetMaxHitsA().
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
48 // processing time in case many noise and/or afterpulse signals are not
49 // recognised by the hit cleaning procedure.
51 // Information about the actual parameter settings can be found in the event
52 // structure itself via the device named "IceDwalk".
54 // The various reconstruction steps are summarised as follows :
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
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().
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.
73 // 2) Each TE will obtain so called associated hits.
74 // A hit is associated to a TE when it fulfills both the conditions
76 // -30 < tres < 300 ns
79 // tres : time residual
80 // Difference between the observed hit time and the time expected
81 // for a direct photon hit.
82 // dhit : Distance traveled by the cherenkov photon from the track to the hit position
83 // lambda : Photon scattering length in ice
85 // By default F is set to 3.07126, but this can be modified via the memberfunction
88 // 3) Construction of track candidates (TC's).
89 // These are TE's that fulfill both the conditions
94 // where we have defined :
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.
100 // The track quality number qtc is defined as follows :
102 // qtc=nah*(term1+term2)-term3-term4-term5
104 // here we have defined :
106 // term1=2*spread/span
107 // term2=2*spreadL/spanL
108 // term3=|spread-expspread|/spread
109 // term4=|spreadL-expspreadL|/spreadL
110 // term5=|medianT|/spreadT
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|.
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| >
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.
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).
136 // 4) The remaining track candidates are clustered into jets when their directions
137 // are within a certain maximum opening angle.
138 // In addition a track candidate must within a certain maximum distance
139 // of the jet starting TC in order to get clustered.
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.
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.
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.
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 :
159 // qtcjet=|jet momentum|/ntracks
161 // This jet quality number is refined on basis of the number of hits
162 // associated to the jet as :
164 // qtcjet=qtcjet+0.2*(nah-nahmax)
166 // where we have defined :
168 // nah : Number of associated hits for the specific jet.
169 // nahmax : Maximum number of associated hits encountered for the jets.
171 // This qtcjet value is then used to order the various jets w.r.t.
172 // decreasing qtcjet quality number.
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.
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.
183 // The latter criterion prevents merging of (nearly) parallel tracks/jets
184 // crossing the detector a very different locations (e.g. muon bundles).
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.
189 // The default maximum opening angle is half the TC maximum opening angle,
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.
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.
200 // Note : Setting the maximum jet opening angle to <=0 will prevent
201 // the merging of jets.
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.
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.
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.
214 // The reference point data of the jet will provide the r0 and t0
215 // (i.e. reference point) of the track.
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.
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).
231 // Note : In case the maximum jet opening angle was specified <0,
232 // only the jet with the highest quality number will appear
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.
237 // For further details the user is referred to NIM A524 (2004) 169.
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).
243 //--- Author: Nick van Eijndhoven 07-oct-2005 Utrecht University
244 //- Modified: NvE $Date$ Utrecht University
245 ///////////////////////////////////////////////////////////////////////////
247 #include "IceDwalk.h"
248 #include "Riostream.h"
250 ClassImp(IceDwalk) // Class implementation to enable ROOT I/O
252 IceDwalk::IceDwalk(const char* name,const char* title) : TTask(name,title)
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.
267 fJangmax=fTangmax/2.;
278 ///////////////////////////////////////////////////////////////////////////
279 IceDwalk::~IceDwalk()
281 // Default destructor.
283 ///////////////////////////////////////////////////////////////////////////
284 void IceDwalk::SetDmin(Float_t d)
286 // Set minimum hit distance (in m) to form a track element.
287 // In the constructor the default has been set to 75 meter.
290 ///////////////////////////////////////////////////////////////////////////
291 void IceDwalk::SetDtmarg(Int_t dt)
293 // Set maximum hit time difference margin (in ns) for track elements.
294 // In the constructor the default has been set to 0 ns.
297 ///////////////////////////////////////////////////////////////////////////
298 void IceDwalk::SetMaxDhit(Float_t d)
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.
304 ///////////////////////////////////////////////////////////////////////////
305 void IceDwalk::SetTangmax(Float_t ang)
307 // Set maximum angular separation (in deg) for track candidate clustering
309 // In the constructor the default has been set to 15 deg, in accordance
310 // to NIM A524 (2004) 180.
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.
320 ///////////////////////////////////////////////////////////////////////////
321 void IceDwalk::SetTdistmax(Float_t d,Int_t invol)
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.
338 ///////////////////////////////////////////////////////////////////////////
339 void IceDwalk::SetJangmax(Float_t ang,Int_t iter)
341 // Set angular separation (in deg) within which jets are merged into 1
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.
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.
367 ///////////////////////////////////////////////////////////////////////////
368 void IceDwalk::SetJdistmax(Float_t d,Int_t invol)
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.
384 ///////////////////////////////////////////////////////////////////////////
385 void IceDwalk::SetMaxModA(Int_t nmax)
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,
394 // which implies no selection on maximum module multiplicity.
395 // See also the memberfunction SetMinModA().
398 ///////////////////////////////////////////////////////////////////////////
399 void IceDwalk::SetMinModA(Int_t nmin)
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().
409 ///////////////////////////////////////////////////////////////////////////
410 void IceDwalk::SetMaxHitsA(Int_t nmax)
412 // Set the maximum number of good hits per Amanda module to be processed.
415 // nmax = 0 : No maximum limit set; all good hits will be processed
416 // < 0 : No hits will be processed
418 // In case the user selects a maximum number of good hits per module, all the
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...
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.
425 ///////////////////////////////////////////////////////////////////////////
426 void IceDwalk::SetVgroupUsage(Int_t flag)
428 // (De)activate the distinction between v_phase and v_group of the Cherenkov light.
430 // flag = 0 : No distinction between v_phase and v_group
431 // = 1 : Separate treatment of v_phase and v_group
433 // By default the distinction between v_phase and v_group is activated
434 // in the constructor of this class.
437 ///////////////////////////////////////////////////////////////////////////
438 void IceDwalk::SetTrackName(TString s)
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.
443 // By default the produced first guess tracks have the name of the class
444 // by which they were produced.
447 ///////////////////////////////////////////////////////////////////////////
448 void IceDwalk::SetCharge(Float_t charge)
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.
456 ///////////////////////////////////////////////////////////////////////////
457 void IceDwalk::Exec(Option_t* opt)
459 // Implementation of the direct walk track reconstruction.
462 AliJob* parent=(AliJob*)(gROOT->GetListOfTasks()->FindObject(name.Data()));
466 fEvt=(IceEvent*)parent->GetObject("IceEvent");
469 // Enter the reco parameters as a device in the event
471 params.SetNameTitle(ClassName(),"Reco parameters");
472 params.SetSlotName("Dmin",1);
473 params.SetSlotName("Dtmarg",2);
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);
487 params.SetSignal(fDmin,1);
488 params.SetSignal(fDtmarg,2);
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);
502 fEvt->AddDevice(params);
504 if (fMaxhitsA<0) return;
506 // Fetch all fired Amanda OMs for this event
507 TObjArray* aoms=fEvt->GetDevices("IceAOM");
509 Int_t naoms=aoms->GetEntries();
512 // Check for the minimum and/or maximum number of good fired Amanda OMs
514 for (Int_t iom=0; iom<naoms; iom++)
516 IceGOM* omx=(IceGOM*)aoms->At(iom);
518 if (omx->GetDeadValue("ADC") || omx->GetDeadValue("LE") || omx->GetDeadValue("TOT")) continue;
521 if (ngood<fMinmodA || ngood>fMaxmodA) return;
523 const Float_t c=0.299792458; // Light speed in vacuum in meters per ns
525 // Storage of track elements.
539 Float_t dist=0,dist2=0;
545 // Check the hits of Amanda OM pairs for possible track elements.
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.
549 for (Int_t i1=0; i1<naoms; i1++) // First OM of the pair
551 IceGOM* omx1=(IceGOM*)aoms->At(i1);
553 if (omx1->GetDeadValue("LE")) continue;
554 r1=omx1->GetPosition();
555 // Select all the good hits of this first OM
557 // Determine the max. number of hits to be processed for this OM
559 if (fMaxhitsA>0 && omx1->GetNhits()>fMaxhitsA) ordered=omx1->SortHits("LE",1,0,7);
561 for (Int_t j1=1; j1<=omx1->GetNhits(); j1++)
565 if (nh1>=fMaxhitsA) break;
566 sx1=(AliSignal*)ordered->At(j1-1);
570 sx1=omx1->GetHit(j1);
573 if (sx1->GetDeadValue("ADC") || sx1->GetDeadValue("LE") || sx1->GetDeadValue("TOT")) continue;
575 // Also store all good hits in the total hit array
580 // No further pair to be formed with the last OM in the list
581 if (i1==(naoms-1)) break;
583 nh1=hits1.GetEntries();
586 for (Int_t i2=i1+1; i2<naoms; i2++) // Second OM of the pair
588 IceGOM* omx2=(IceGOM*)aoms->At(i2);
590 if (omx2->GetDeadValue("LE")) continue;
591 r2=omx2->GetPosition();
595 if (dist<=fDmin) continue;
597 // Select all the good hits of this second OM
599 // Determine the max. number of hits to be processed for this OM
601 if (fMaxhitsA>0 && omx2->GetNhits()>fMaxhitsA) ordered=omx2->SortHits("LE",1,0,7);
603 for (Int_t j2=1; j2<=omx2->GetNhits(); j2++)
607 if (nh2>=fMaxhitsA) break;
608 sx2=(AliSignal*)ordered->At(j2-1);
612 sx2=omx2->GetHit(j2);
615 if (sx2->GetDeadValue("ADC") || sx2->GetDeadValue("LE") || sx2->GetDeadValue("TOT")) continue;
620 nh2=hits2.GetEntries();
623 // Position r0 in between the two OMs and normalised relative direction r12
625 r0.SetPosition((Ali3Vector&)rsum);
628 // Check all hit pair combinations of these two OMs for possible track elements
629 dtmax=dist/c+float(fDtmarg);
630 for (Int_t ih1=0; ih1<nh1; ih1++) // Hits of first OM
632 sx1=(AliSignal*)hits1.At(ih1);
634 for (Int_t ih2=0; ih2<nh2; ih2++) // Hits of second OM
636 sx2=(AliSignal*)hits2.At(ih2);
638 t1=sx1->GetSignal("LE",7);
639 t2=sx2->GetSignal("LE",7);
643 if (fabs(dt)>=dtmax) continue;
648 r0.SetTimestamp((AliTimestamp&)*fEvt);
649 AliTimestamp* tsx=r0.GetTimestamp();
650 tsx->Add(0,0,(int)t0);
651 te->SetReferencePoint(r0);
652 te->Set3Momentum(r12);
655 } // end of loop over the second OM of the pair
656 } // end of loop over first OM of the pair
658 // Association of hits to the various track elements.
661 AssociateHits(tes,hits,qmax,nahmax);
663 // Selection on quality (Q value) in case of multiple track candidates
664 SelectQvalue(tes,qmax);
666 Int_t nte=tes.GetEntries();
669 // Clustering of track candidates into jets
672 ClusterTracks(tes,jets,qmax);
674 Int_t njets=jets.GetEntries();
677 // Order the jets w.r.t. decreasing quality value
678 ordered=fEvt->SortJets(-2,&jets);
679 TObjArray jets2(*ordered);
684 // Production and storage of the final tracks
687 ///////////////////////////////////////////////////////////////////////////
688 void IceDwalk::AssociateHits(TObjArray& tes,TObjArray& hits,Float_t& qmax,Int_t& nahmax)
690 // Association of hits to the various track elements.
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)
699 // Angular reduction of complement of thetac due to v_phase and v_group difference
701 if (fVgroup) alphac=atan((1.-npice/ngice)/sqrt(npice*npice-1.));
703 Int_t nte=tes.GetEntries();
704 Int_t nh=hits.GetEntries();
711 Float_t dist,t0,tgeo,tres;
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");
742 Int_t nah; // Number of associated hits for a certain TE
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;
749 for (Int_t jte=0; jte<nte; jte++)
751 AliTrack* te=(AliTrack*)tes.At(jte);
753 AliPosition* tr0=te->GetReferencePoint();
754 AliTimestamp* tt0=tr0->GetTimestamp();
755 t0=fEvt->GetDifference(tt0,"ns");
756 p=te->Get3Momentum();
760 for (Int_t jh=0; jh<nh; jh++)
762 AliSignal* sx1=(AliSignal*)hits.At(jh);
764 IceGOM* omx=(IceGOM*)sx1->GetDevice();
766 r1=omx->GetPosition();
767 d=te->GetDistance(r1);
770 dist=hproj+d/tan(pi/2.-thetac-alphac);
772 t1=sx1->GetSignal("LE",7);
775 d=d/sin(thetac); // The distance traveled by a cherenkov photon
777 if (tres<-30 || tres>300 || d>fMaxdhit*lambda) continue;
779 // Associate this hit to the TE
781 levers.Enter(fabs(hproj));
786 // Determine the Q quality of the various TE's.
787 // Good quality TE's will be called track candidates (TC's)
788 nah=te->GetNsignals();
789 if (nah>nahmax) nahmax=nah;
790 lmin=levers.GetMinimum(1);
791 lmax=levers.GetMaximum(1);
793 medianl=levers.GetMedian(1);
794 meanl=levers.GetMean(1);
795 sigmal=levers.GetSigma(1);
796 spreadl=levers.GetSpread(1);
797 // Expected spread for a flat distribution
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
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);
816 if (span>0) term1=2.*spread/span;
819 if (spanl>0) term2=2.*spreadl/spanl;
822 if (spread>0) term3=fabs(spread-expspread)/spread;
825 if (spreadl>0) term4=fabs(spreadl-expspreadl)/spreadl;
828 if (spreadt>0) term5=fabs(mediant)/spreadt;
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
833 fit.SetSignal(qtc,"QTC");
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");
855 te->SetFitDetails(fit);
856 if (qtc>qmax) qmax=qtc;
859 ///////////////////////////////////////////////////////////////////////////
860 void IceDwalk::SelectQvalue(TObjArray& tes,Float_t qmax)
862 // Perform selection on Q value in case of multiple track candidates
864 Int_t nte=tes.GetEntries();
868 for (Int_t jtc=0; jtc<nte; jtc++)
870 AliTrack* te=(AliTrack*)tes.At(jtc);
872 nah=te->GetNsignals();
873 AliSignal* sx1=(AliSignal*)te->GetFitDetails();
875 if (sx1) qtc=sx1->GetSignal("QTC");
876 if (!nah || qtc<0.8*qmax)
881 else // Set Q value as momentum to provide a weight for jet clustering
885 p=te->Get3Momentum();
893 ///////////////////////////////////////////////////////////////////////////
894 void IceDwalk::ClusterTracks(TObjArray& tes,TObjArray& jets,Float_t qmax)
896 // Cluster track candidates within a certain opening angle into jets.
897 // Also the track should be within a certain maximum distance of the
898 // starting track in order to get clustered.
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.
904 Int_t nte=tes.GetEntries();
908 Float_t vec[3],err[3];
910 Float_t t0,dist,dist2;
911 Int_t nah=0,nahmax=0; // Determine the max. number of associated hits for the jets
913 for (Int_t jtc1=0; jtc1<nte; jtc1++)
915 AliTrack* te=(AliTrack*)tes.At(jtc1);
917 AliPosition* x1=te->GetReferencePoint();
919 AliTimestamp* ts1=x1->GetTimestamp();
921 AliJet* jx=new AliJet();
925 x1->GetPosition(vec,"car");
926 pos.Enter(vec[0],vec[1],vec[2]);
927 t0=fEvt->GetDifference(ts1,"ns");
929 for (Int_t jtc2=0; jtc2<nte; jtc2++)
931 if (jtc2==jtc1) continue;
932 AliTrack* te2=(AliTrack*)tes.At(jtc2);
934 ang=te->GetOpeningAngle(*te2,"deg");
937 AliPosition* x2=te2->GetReferencePoint();
939 AliTimestamp* ts2=x2->GetTimestamp();
943 dist=te->GetDistance(te2);
947 dist=te->GetDistance(x2);
948 dist2=te2->GetDistance(x1);
949 if (dist2<dist) dist=dist2;
953 x2->GetPosition(vec,"car");
954 pos.Enter(vec[0],vec[1],vec[2]);
955 t0=fEvt->GetDifference(ts2,"ns");
962 // Set the reference point data for this jet
963 for (Int_t j=1; j<=3; j++)
965 vec[j-1]=pos.GetMean(j);
966 err[j-1]=pos.GetSigma(j);
968 r0.SetPosition(vec,"car");
969 r0.SetPositionErrors(err,"car");
970 r0.SetTimestamp((AliTimestamp&)*fEvt);
971 AliTimestamp* jt0=r0.GetTimestamp();
973 jt0->Add(0,0,(int)t0);
974 jx->SetReferencePoint(r0);
976 // Store this jet for further processing if ntracks>1
977 if (jx->GetNtracks() > 1 || fTangmax<=0)
980 nah=jx->GetNsignals();
981 if (nah>nahmax) nahmax=nah;
983 else // Only keep single-track jets which have qtc=qmax
985 AliSignal* sx1=(AliSignal*)te->GetFitDetails();
987 if (sx1) qtc=sx1->GetSignal("QTC");
988 if (qtc>=(qmax-1.e-10))
991 nah=jx->GetNsignals();
992 if (nah>nahmax) nahmax=nah;
1001 Int_t njets=jets.GetEntries();
1004 // The sum of 0.15*(nah-nahmax) and average qtc value per track for each jet
1005 // will be stored as the jet energy to enable sorting on this value lateron
1008 for (Int_t ijet=0; ijet<njets; ijet++)
1010 AliJet* jx=(AliJet*)jets.At(ijet);
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);
1019 ///////////////////////////////////////////////////////////////////////////
1020 void IceDwalk::MergeJets(TObjArray& jets2)
1022 // Merge jets within a certain opening to provide the final track(s).
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.
1030 Int_t njets=jets2.GetEntries();
1034 Int_t ntk,nah,nahmax;
1035 Float_t ang,dist,dist2,t0;
1039 Float_t vec[3],err[3];
1047 for (Int_t jet1=0; jet1<njets; jet1++)
1049 jx1=(AliJet*)jets2.At(jet1);
1051 AliPosition* x1=jx1->GetReferencePoint();
1053 AliTimestamp* ts1=x1->GetTimestamp();
1057 x1->GetPosition(vec,"car");
1058 pos.Enter(vec[0],vec[1],vec[2]);
1059 t0=fEvt->GetDifference(ts1,"ns");
1061 for (Int_t jet2=0; jet2<njets; jet2++)
1063 jx2=(AliJet*)jets2.At(jet2);
1064 if (!jx2 || jet2==jet1) continue;
1065 AliPosition* x2=jx2->GetReferencePoint();
1067 AliTimestamp* ts2=x2->GetTimestamp();
1069 ang=jx1->GetOpeningAngle(*jx2,"deg");
1074 dist=jx1->GetDistance(jx2);
1078 dist=jx1->GetDistance(x2);
1079 dist2=jx2->GetDistance(x1);
1080 if (dist2<dist) dist=dist2;
1082 if (dist<=fJdistmax)
1084 x2->GetPosition(vec,"car");
1085 pos.Enter(vec[0],vec[1],vec[2]);
1086 t0=fEvt->GetDifference(ts2,"ns");
1088 for (Int_t jtk=1; jtk<=jx2->GetNtracks(); jtk++)
1090 AliTrack* te=jx2->GetTrack(jtk);
1094 jets2.RemoveAt(jet2);
1095 if (fJiterate) merged=1;
1098 } // End of jet2 loop
1100 // Set the reference point data for this jet
1101 for (Int_t k=1; k<=3; k++)
1103 vec[k-1]=pos.GetMean(k);
1104 err[k-1]=pos.GetSigma(k);
1106 r0.SetPosition(vec,"car");
1107 r0.SetPositionErrors(err,"car");
1108 r0.SetTimestamp((AliTimestamp&)*fEvt);
1109 AliTimestamp* jt0=r0.GetTimestamp();
1111 jt0->Add(0,0,(int)t0);
1112 jx1->SetReferencePoint(r0);
1114 nah=jx1->GetNsignals();
1115 if (nah>nahmax) nahmax=nah;
1116 } // End of jet1 loop
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
1122 for (Int_t jjet=0; jjet<njets; jjet++)
1124 AliJet* jx=(AliJet*)jets2.At(jjet);
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);
1133 // Order the jets w.r.t. decreasing quality value
1134 TObjArray* ordered=fEvt->SortJets(-2,&jets2);
1135 njets=ordered->GetEntries();
1137 for (Int_t icopy=0; icopy<njets; icopy++)
1139 jets2.Add(ordered->At(icopy));
1141 } // End of iterative while loop
1144 ///////////////////////////////////////////////////////////////////////////
1145 void IceDwalk::StoreTracks(TObjArray& jets2)
1147 // Store every jet as a reconstructed track in the event structure.
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.
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.
1158 Int_t njets=jets2.GetEntries();
1160 if (fTrackname=="") fTrackname=ClassName();
1161 TString title=ClassName();
1162 title+=" reco track";
1164 t.SetNameTitle(fTrackname.Data(),title.Data());
1165 t.SetCharge(fCharge);
1167 for (Int_t jet=0; jet<njets; jet++)
1169 AliJet* jx=(AliJet*)jets2.At(jet);
1171 AliPosition* ref=jx->GetReferencePoint();
1174 AliTrack* trk=fEvt->GetTrack(fEvt->GetNtracks());
1176 trk->SetId(fEvt->GetNtracks(1)+1);
1177 p=jx->Get3Momentum();
1179 trk->Set3Momentum(p);
1180 trk->SetReferencePoint(*ref);
1181 AliTimestamp* tt0=ref->GetTimestamp();
1182 if (tt0) trk->SetTimestamp(*tt0);
1183 for (Int_t jt=1; jt<=jx->GetNtracks(); jt++)
1185 AliTrack* tx=jx->GetTrack(jt);
1187 for (Int_t is=1; is<=tx->GetNsignals(); is++)
1189 AliSignal* sx1=tx->GetSignal(is);
1190 if (sx1) sx1->AddTrack(*trk);
1194 // Only take the jet with the highest quality number
1195 // (i.e. the first jet in the list) when the user had selected
1196 // this reconstruction mode.
1197 if (fJangmax<0) break;
1200 ///////////////////////////////////////////////////////////////////////////