warning fixed
[u/mrichter/AliRoot.git] / RALICE / AliHelix.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 // $Id$
17
18 ///////////////////////////////////////////////////////////////////////////
19 // Class AliHelix
20 // Representation and extrapolation of AliTracks in a magnetic field.
21 //
22 // This class is meant to provide a means to display and extrapolate
23 // AliTrack objects in the presence of a constant homogeneous magnetic field. 
24 //
25 // For track/event displays the line width, colour etc... can be set using the
26 // standard facilities (see TAttLine).
27 // By default the linewith is set to 2 and the colour set to -1 in the constructor.
28 // The latter results in an automatic colour coding according to the track charge
29 // with the convention positive=red neutral=green negative=blue.
30 //
31 // To indicate the track starting point, the memberfunction SetMarker()
32 // may be used.
33 // By default no marker will be displayed. 
34 //
35 // Examples :
36 // ==========
37 //
38 // Display and extrapolation of individual tracks 
39 // ----------------------------------------------
40 // Float_t vec[3];
41 // AliPosition r1;
42 // Ali3Vector p;
43 // AliTrack t;
44 //
45 // vec[0]=0;
46 // vec[1]=0;
47 // vec[2]=0;
48 // r1.SetVector(vec,"car");
49 //
50 // vec[0]=1;
51 // vec[1]=0;
52 // vec[2]=0.3;
53 // p.SetVector(vec,"car");
54 //
55 // t.Set3Momentum(p);
56 // t.SetBeginPoint(r1);
57 // t.SetCharge(-1);
58 // t.SetMass(0.139);
59 //
60 // // The magnetic field vector in Tesla
61 // Ali3Vector b;
62 // vec[0]=0;
63 // vec[1]=0;
64 // vec[2]=1;
65 // b.SetVector(vec,"car");
66 //
67 // AliHelix* helix=new AliHelix();
68 // helix->SetB(b);
69 // helix->SetTofmax(1e-7);
70 //
71 // TCanvas* c1=new TCanvas("c1","c1");
72 // TView* view=new TView(1);
73 // view->SetRange(-1000,-1000,-1000,1000,1000,1000);
74 // view->ShowAxis();
75 //
76 // // Track displays 
77 // Double_t range[2]={0,600};
78 // helix->Display(&t,range,3);
79 // t.SetCharge(-t.GetCharge());
80 // helix->Display(&t);
81 //
82 // // Track extrapolation
83 // Double_t pars[3]={550,0.001,3};
84 // AliPosition* rext=helix->Extrapolate(&t,pars);
85 // if (rext) rext->Data();
86 // ======================================================================
87 //
88 // Online display of events generated via AliCollider
89 // -------------------------------------------------- 
90 // Int_t nevents=5;   // Number of events to be generated
91 // Int_t jrun=1;      // The run number of this batch of generated events
92 //
93 // cout << " ***" << endl;
94 // cout << " *** AliCollider run for " << nevents << " events." << endl; 
95 // cout << " ***" << endl;
96 //
97 // AliCollider* gen=new AliCollider();
98 //
99 // gen->OpenFortranFile(6,"dump.log");
100 //
101 // gen->SetVertexMode(2);
102 // gen->SetResolution(1e-4);
103 //
104 // gen->SetRunNumber(jrun);
105 // gen->SetPrintFreq(1);
106 //
107 // gen->SetSpectatorPmin(0.01);
108 //
109 // Int_t zp=1;
110 // Int_t ap=1;
111 // Int_t zt=2;
112 // Int_t at=4;
113 //
114 // gen->Init("fixt",zp,ap,zt,at,158);
115 //
116 // AliHelix* helix=new AliHelix();
117 // Float_t vec[3]={0,2,0};
118 // Ali3Vector b;
119 // b.SetVector(vec,"car");
120 // helix->SetB(b);
121 //
122 // helix->Refresh(-1); // Refresh display after each event
123 //
124 // TCanvas* c1=new TCanvas("c1","c1");
125 // TView* view=new TView(1);
126 // view->SetRange(-200,-200,-200,200,200,200);
127 // view->ShowAxis();
128 //
129 // // Prepare random number sequence for this run
130 // // to obtain the number of participants for each event
131 // AliRandom rndm(abs(jrun));
132 // Float_t* rans=new Float_t[nevents];
133 // rndm.Uniform(rans,nevents,2,ap+at);
134 // Int_t npart=0;
135 // Int_t ntk=0;
136 // for (Int_t i=0; i<nevents; i++)
137 // {
138 //  npart=rans[i];
139 //  gen->MakeEvent(npart);
140 //  AliEvent* evt=gen->GetEvent();
141 //  if (evt)
142 //  {
143 //   helix->Display(evt);
144 //   c1->Update();
145 //   gSystem->Sleep(5000); // Some delay to keep the display on screen
146 //  }
147 // }
148 // ======================================================================
149 //
150 //--- Author: Nick van Eijndhoven 17-jun-2004 Utrecht University
151 //- Modified: NvE $Date$ Utrecht University
152 ///////////////////////////////////////////////////////////////////////////
153
154 #include <cstdlib>
155 #include "AliHelix.h"
156 #include "Riostream.h"
157  
158 ClassImp(AliHelix) // Class implementation to enable ROOT I/O
159  
160 AliHelix::AliHelix() : THelix()
161 {
162 // Default constructor
163  fRefresh=0;
164  fCurves=0;
165  fExt=0;
166  fTofmax=1e-8;
167  fMstyle=-1;
168  fMsize=0;
169  fMcol=0;
170  fEnduse=1;
171
172  fLineWidth=2;
173  fLineColor=-1;
174 }
175 ///////////////////////////////////////////////////////////////////////////
176 AliHelix::~AliHelix()
177 {
178 // Destructor to delete dynamically allocated memory.
179  if (fCurves)
180  {
181   delete fCurves;
182   fCurves=0;
183  }
184  if (fExt)
185  {
186   delete fExt;
187   fExt=0;
188  }
189 }
190 ///////////////////////////////////////////////////////////////////////////
191 AliHelix::AliHelix(const AliHelix& h) : THelix(h)
192 {
193 // Copy constructor
194  fB=h.fB;
195  fRefresh=h.fRefresh;
196  fTofmax=h.fTofmax;
197  fMstyle=h.fMstyle;
198  fMsize=h.fMsize;
199  fMcol=h.fMcol;
200  fEnduse=h.fEnduse;
201 }
202 ///////////////////////////////////////////////////////////////////////////
203 void AliHelix::SetB(Ali3Vector& b)
204 {
205 // Set the magnetic field vector in Tesla.
206  fB=b;
207
208  if (fB.GetNorm()>0)
209  {
210   Double_t axis[3];
211   fB.GetVector(axis,"car");
212   SetAxis(axis);
213  }
214 }
215 ///////////////////////////////////////////////////////////////////////////
216 Ali3Vector& AliHelix::GetB()
217 {
218 // Provide the magnetic field vector in Tesla.
219  return fB;
220 }
221 ///////////////////////////////////////////////////////////////////////////
222 void AliHelix::SetTofmax(Float_t tof)
223 {
224 // Set the maximum time of flight for straight tracks in seconds.
225 // This maximum tof will be used for drawing etc... in case no begin
226 // and endpoints can be determined from the track info.
227 // Notes :
228 // -------
229 // 1) In case the user specifies an explicit range, it will override
230 //    the maximum tof limit.
231 // 2) By default the tofmax is set to 10 ns in the AliHelix constructor.
232  fTofmax=tof;
233 }
234 ///////////////////////////////////////////////////////////////////////////
235 Float_t AliHelix::GetTofmax() const
236 {
237 // Provide the maximum time of flight for straight tracks in seconds.
238  return fTofmax;
239 }
240 ///////////////////////////////////////////////////////////////////////////
241 void AliHelix::SetMarker(Int_t style,Float_t size,Int_t col)
242 {
243 // Specify the marker (style, size and colour) to indicate the starting point
244 // of a track in a display.
245 // In case col<0 the marker will have the same color as the track itself.
246 // 
247 // Defaults are style=8, size=0.2 and col=-1.
248  
249  fMstyle=style;
250  fMsize=size;
251  fMcol=col;
252 }
253 ///////////////////////////////////////////////////////////////////////////
254 void AliHelix::UseEndPoint(Int_t mode)
255 {
256 // Select usage of track endpoint in drawing and extrapolation.
257 // This allows correct event displays even for very long tracks.
258 //
259 // mode = 0 : Do not use the track endpoint
260 //        1 : Use the track endpoint
261 // 
262 // The default value is mode=1 (which is also set in the constructor).
263
264  if (mode==0 || mode==1) fEnduse=mode; 
265 }
266 ///////////////////////////////////////////////////////////////////////////
267 void AliHelix::MakeCurve(AliTrack* t,Double_t* range,Int_t iaxis,Double_t scale)
268 {
269 // Make the helix curve for the specified AliTrack.
270 // Detailed information of all the helix points can be obtained via the
271 // GetN() and GetP() memberfunctions of TPolyLine3D.
272 // In case one wants to display or extrapolate an AliTrack it is preferable
273 // to use the Display() or Extrapolate() memberfunctions.
274 // It is assumed that the track charge is stored in elementary units
275 // (i.e. charge=1 for a proton).
276 // The input argument "scale" specifies the unit scale for the various
277 // locations where scale=0.01 indicates unit scales in cm etc...
278 // In case scale<=0, the unit scale for locations is determined from the
279 // begin, reference or endpoint of the track. If neither of these
280 // positions is present, all locations are assumed to be given in meter.
281 // The lower and upper bounds for the range are specified by range[0] and
282 // range[1] and the argument "iaxis" indicates along which axis this range
283 // is specified.
284 // The range can be specified either in the LAB frame or in the Helix frame.
285 // The latter is the frame in which the Z axis points in the B direction.
286 //
287 // The conventions for the "iaxis" argument are the following :
288 // iaxis = 1 ==> X axis in the LAB frame
289 //         2 ==> Y axis in the LAB frame
290 //         3 ==> Z axis in the LAB frame
291 //        -1 ==> X axis in the Helix frame
292 //        -2 ==> Y axis in the Helix frame
293 //        -3 ==> Z axis in the Helix frame
294 //
295 // In case range=0 the begin/end/reference points of the AliTrack and the
296 // maximum time of flight (see the SetTofmax() memberfunction) will be used
297 // and an appropriate choice for the iaxis parameter will be made automatically
298 // based on the track kinematics.
299 // In case the reference point is not present, the begin or endpoint will be used
300 // as reference point for the 3-momentum specification. If neither of these positions
301 // is present, (0,0,0) will be taken as the reference point.
302 // 
303 // The default values are range=0, iaxis=3 and scale=-1.
304
305  SetPolyLine(0); // Reset the polyline data points
306
307  if (!t || (range && !iaxis)) return;
308
309  Double_t energy=t->GetEnergy(1); // Track energy in GeV
310  Double_t betanorm=t->GetBeta();
311
312  if (energy<=0 || betanorm<=0) return;
313
314  AliPosition* rbeg=t->GetBeginPoint();
315  AliPosition* rend=0;
316  if (fEnduse) rend=t->GetEndPoint();
317  AliPosition* rref=t->GetReferencePoint();
318
319  // Magnetic field vector or default Z-direction
320  Double_t bvec[3]={0,0,1};
321  if (fB.GetNorm()>0) fB.GetVector(bvec,"car");
322
323  // The unit scale for locations if not specified by the user
324  if (scale<=0)
325  {
326   scale=1; // Set default to meter
327   if (rbeg)
328   {
329    scale=rbeg->GetUnitScale();
330   }
331   else if (rend)
332   {
333    scale=rend->GetUnitScale();
334   }
335   else if (rref)
336   {
337    scale=rref->GetUnitScale();
338   }
339  }
340
341  Double_t c=2.99792458e8/scale; // Lightspeed in the selected unit scale
342
343  // The helix angular frequency
344  Double_t w=9e7*(t->GetCharge()*fB.GetNorm())/energy;
345
346  // The particle velocity in the LAB frame
347  Ali3Vector beta=t->GetBetaVector();
348  Ali3Vector v=beta*c;
349  Double_t vel[3];
350  v.GetVector(vel,"car");
351
352  // The particle velocity in the Helix frame
353  Ali3Vector betaprim=beta.GetPrimed(fRotMat);
354  v=v.GetPrimed(fRotMat);
355  Double_t velprim[3];
356  v.GetVector(velprim,"car");
357
358  // Check compatibility of velocity and range specification.
359  if (range)
360  {
361   Double_t betavec[3];
362   if (iaxis>0) beta.GetVector(betavec,"car");
363   if (iaxis<0) betaprim.GetVector(betavec,"car");
364   if (fabs(betavec[abs(iaxis)-1])/betanorm<1e-10) return;
365  }
366
367  // The LAB location in which the velocity of the particle is defined
368  Double_t loc[3]={0,0,0};
369  Ali3Vector rx;
370  Double_t scalex=-1;
371  if (rref)
372  {
373   rx=(Ali3Vector)(*rref);
374   scalex=rref->GetUnitScale();
375  }
376  else if (rbeg)
377  {
378   rx=(Ali3Vector)(*rbeg);
379   scalex=rbeg->GetUnitScale();
380  }
381  else if (rend)
382  {
383   rx=(Ali3Vector)(*rend);
384   scalex=rend->GetUnitScale();
385  }
386
387  if (scalex>0 && (scalex/scale>1.1 || scale/scalex>1.1)) rx*=scalex/scale;
388  rx.GetVector(loc,"car");
389
390  // Initialisation of Helix kinematics
391  SetHelix(loc,vel,w,0,kUnchanged,bvec);
392
393  Int_t bend=0;
394  if (fabs(w)>0 && fabs(fVt)>0) bend=1;
395
396  // Flight time boundaries.
397  // The time origin t=0 is chosen to indicate the position in which
398  // the particle velocity was defined.
399  // The total flight time is initialised to the (user specified) tofmax.
400  Double_t tmin=0,tmax=0;
401  Double_t tof=fTofmax;
402  Double_t dum=0;
403
404  // The trajectory begin and end points
405  Double_t vec1[3]={0,0,0};
406  Double_t vec2[3]={0,0,0};
407  Ali3Vector r1;
408  Ali3Vector r2;
409  Double_t scale1=1;
410  Double_t scale2=1;
411
412  if (!bend)
413  {
414   ////////////////////////////////////////
415   // Treatment of straight trajectories //
416   ////////////////////////////////////////
417   Ali3Vector r;
418   if (range) // Specified range allows for exact flight time boundaries
419   {
420    if (iaxis>0)
421    {
422     tmin=(range[0]-loc[iaxis-1])/vel[iaxis-1];
423     tmax=(range[1]-loc[iaxis-1])/vel[iaxis-1];
424    }
425    else
426    {
427     loc[0]=fX0;
428     loc[1]=fY0;
429     loc[2]=fZ0;
430     tmin=(range[0]-loc[abs(iaxis)-1])/velprim[abs(iaxis)-1];
431     tmax=(range[1]-loc[abs(iaxis)-1])/velprim[abs(iaxis)-1];
432    }
433    if (tmax<tmin)
434    {
435     dum=tmin;
436     tmin=tmax;
437     tmax=dum;
438    }
439    // Make the 'curve' in the LAB frame and exit.
440    // Use the parametrisation : r(t)=r0+t*v
441    // using the range based flight time boundaries.
442    // An additional point in the middle of the trajectory is
443    // generated in view of accuracy in the case of extrapolations.
444    tof=tmax-tmin;
445    v=beta*c;
446    r1=rx;
447    r=v*tmin;
448    r1=r1+r;
449    r1.GetVector(vec1,"car");
450    SetNextPoint(float(vec1[0]),float(vec1[1]),float(vec1[2]));
451    r=v*(tof/2.);
452    r2=r1+r;
453    r2.GetVector(vec2,"car");
454    SetNextPoint(float(vec2[0]),float(vec2[1]),float(vec2[2]));
455    r=v*tof;
456    r2=r1+r;
457    r2.GetVector(vec2,"car");
458    SetNextPoint(float(vec2[0]),float(vec2[1]),float(vec2[2]));
459   }
460   else // Automatic range determination
461   {
462    // Initially the point with Z=0 in the Helix frame is taken as a starting point.
463    // In case this point can't be reached, the point in which the particle velocity
464    // was defined is taken as the starting point.
465    // The endpoint is initially obtained by applying the tofmax from the start point.
466    tmin=0;
467    if (fabs(fVz)>0) tmin=-fZ0/fVz;
468    v=beta*c;
469    r1=rx;
470    r=v*tmin;
471    r1=r1+r;
472
473    // Override the initial begin and endpoint settings by the track data
474    if (rbeg)
475    {
476     r1=(Ali3Vector)(*rbeg); 
477     scale1=rbeg->GetUnitScale();
478     // All coordinates in the selected unit scale
479     if (scale1/scale>1.1 || scale/scale1>1.1) r1*=scale1/scale;
480    }
481
482    r=v*fTofmax;
483    r2=r1+r;
484    if (rend)
485    {
486     r2=(Ali3Vector)(*rend); 
487     scale2=rend->GetUnitScale();
488     // All coordinates in the selected unit scale
489     if (scale2/scale>1.1 || scale/scale2>1.1) r2*=scale2/scale;
490    }
491    
492    r1.GetVector(vec1,"car");
493    r2.GetVector(vec2,"car");
494
495    // Make the 'curve' in the LAB frame and exit.
496    SetNextPoint(float(vec1[0]),float(vec1[1]),float(vec1[2]));
497    SetNextPoint(float(vec2[0]),float(vec2[1]),float(vec2[2]));
498   }
499  }
500  else
501  {
502   //////////////////////////////////////
503   // Treatment of curved trajectories //
504   //////////////////////////////////////
505
506   // Initialisation of the flight time boundaries.
507   // Based on the constant motion of the particle along the Helix Z-axis,
508   // the parametrisation z(t)=z0+fVz*t in the Helix frame is used. 
509   // If possible the point with Z=0 in the Helix frame is taken as a starting point.
510   // In case this point can't be reached, the point in which the particle velocity
511   // was defined is taken as the starting point.
512   tmin=0;
513   if (fabs(fVz)>0) tmin=-fZ0/fVz;
514   tmax=tmin+fTofmax;
515
516   if (tmax<tmin)
517   {
518    dum=tmin;
519    tmin=tmax;
520    tmax=dum;
521   }
522
523   // Determination of the range in the helix frame
524
525   if (!range) // Automatic range determination
526   {
527    scale1=1;
528    scale2=1;
529    if (rbeg)
530    {
531     r1=rbeg->GetPrimed(fRotMat);
532     scale1=rbeg->GetUnitScale();
533     // All coordinates in the selected unit scale
534     if (scale1/scale>1.1 || scale/scale1>1.1) r1*=scale1/scale;
535     // Re-calculate the tmin for this new starting point
536     r1.GetVector(vec1,"car");
537     if (fabs(fVz)>0) tmin=(vec1[2]-fZ0)/fVz;
538     tmax=tmin+fTofmax;
539    }
540    if (rend)
541    {
542     r2=rend->GetPrimed(fRotMat);
543     scale2=rend->GetUnitScale();
544     // All coordinates in the selected unit scale
545     if (scale2/scale>1.1 || scale/scale2>1.1) r2*=scale2/scale;
546     r2.GetVector(vec2,"car");
547     if (fabs(fVz)>0) tmax=(vec2[2]-fZ0)/fVz;
548    }
549    // Make the curve on basis of the flight time boundaries and exit
550    if (tmax<tmin)
551    {
552     dum=tmin;
553     tmin=tmax;
554     tmax=dum;
555    }
556    SetRange(tmin,tmax,kHelixT);
557   }
558   else // User explicitly specified range
559   {
560    vec1[abs(iaxis)-1]=range[0];
561    vec2[abs(iaxis)-1]=range[1];
562    r1.SetVector(vec1,"car");
563    r2.SetVector(vec2,"car");
564    if (iaxis>0) // Range specified in LAB frame
565    {
566     r1=r1.GetPrimed(fRotMat);
567     r1.GetVector(vec1,"car");
568     r2=r2.GetPrimed(fRotMat);
569     r2.GetVector(vec2,"car");
570    } 
571    // Determination of the axis component with the
572    // largest range difference
573    Double_t dmax=0;
574    Int_t imax=0;
575    Double_t test=0;
576    for (Int_t i=0; i<3; i++)
577    {
578     test=fabs(vec1[i]-vec2[i]);
579     if (test>dmax)
580     {
581      dmax=test;
582      imax=i;
583     }
584    }
585
586    Double_t rmin=vec1[imax];
587    Double_t rmax=vec2[imax];
588    if (rmax<rmin)
589    {
590     dum=rmin;
591     rmin=rmax;
592     rmax=dum;
593    }
594
595    // The kinematic range boundaries in the helix frame
596    Double_t xmin=fX0-fVt/fW;
597    Double_t xmax=fX0+fVt/fW;
598    Double_t ymin=fY0-fVt/fW;
599    Double_t ymax=fY0+fVt/fW;
600
601    if (xmax<xmin)
602    {
603     dum=xmin;
604     xmin=xmax;
605     xmax=dum;
606    }
607    if (ymax<ymin)
608    {
609     dum=ymin;
610     ymin=ymax;
611     ymax=dum;
612    }
613
614    // Set the range for the helix
615    if (imax==2 && dmax>0) SetRange(rmin,rmax,kHelixZ);
616    if (imax==1)
617    {
618     // Limit range to kinematic boundaries if needed
619     if (rmin<=ymin) rmin=ymin+1e-6*dmax;
620     if (rmax>=ymax) rmax=ymax-1e-6*dmax;
621     if (rmin<rmax) SetRange(rmin,rmax,kHelixY);
622    }
623    if (imax==0)
624    {
625     // Limit range to kinematic boundaries if needed
626     if (rmin<=xmin) rmin=xmin+1e-6*dmax;
627     if (rmax>=xmax) rmax=xmax-1e-6*dmax;
628     if (rmin<rmax) SetRange(rmin,rmax,kHelixX);
629    }
630   }
631  }
632  return;
633 }
634 ///////////////////////////////////////////////////////////////////////////
635 void AliHelix::Display(AliTrack* t,Double_t* range,Int_t iaxis,Double_t scale)
636 {
637 // Display the helix curve of an AliTrack.
638 // Various curves can be displayed together or individually; please refer to
639 // the memberfunction Refresh() for further details.
640 // It is assumed that the track charge is stored in elementary units
641 // (i.e. charge=1 for a proton).
642 // The input argument "scale" specifies the unit scale for the various
643 // locations where scale=0.01 indicates unit scales in cm etc...
644 // In case scale<=0, the unit scale for locations is determined from the
645 // begin, reference or endpoint of the track. If neither of these
646 // positions is present, all locations are assumed to be given in meter.
647 // The lower and upper bounds for the range are specified by range[0] and
648 // range[1] and the argument "iaxis" indicates along which axis this range
649 // is specified.
650 // The range can be specified either in the LAB frame or in the Helix frame.
651 // The latter is the frame in which the Z axis points in the B direction.
652 //
653 // The conventions for the "iaxis" argument are the following :
654 // iaxis = 1 ==> X axis in the LAB frame
655 //         2 ==> Y axis in the LAB frame
656 //         3 ==> Z axis in the LAB frame
657 //        -1 ==> X axis in the Helix frame
658 //        -2 ==> Y axis in the Helix frame
659 //        -3 ==> Z axis in the Helix frame
660 //
661 // In case range=0 the begin/end/reference points of the AliTrack and the
662 // maximum time of flight (see the SetTofmax() memberfunction) will be used
663 // and an appropriate choice for the iaxis parameter will be made automatically
664 // based on the track kinematics.
665 // In case the reference point is not present, the begin or endpoint will be used
666 // as reference point for the 3-momentum specification. If neither of these positions
667 // is present, (0,0,0) will be taken as the reference point.
668 // 
669 // The default values are range=0, iaxis=3 and scale=-1.
670 //
671 // Note :
672 // ------
673 // Before any display activity, a TCanvas and a TView have to be initiated
674 // first by the user like for instance
675 // 
676 // TCanvas* c1=new TCanvas("c1","c1");
677 // TView* view=new TView(1);
678 // view->SetRange(-1000,-1000,-1000,1000,1000,1000);
679 // view->ShowAxis();
680 //
681 // The user can also use the 3D viewing facilities from the TCanvas menu
682 // to open an appropriate view. 
683
684  if (!t || (range && !iaxis)) return;
685
686  MakeCurve(t,range,iaxis,scale);
687
688  if (fRefresh>0) Refresh(fRefresh);
689
690  Int_t np=GetLastPoint()+1;
691  if (!np) return;
692
693  Float_t* points=GetP();
694  TPolyLine3D* curve=new TPolyLine3D(np,points);
695
696  curve->SetLineWidth(fLineWidth);
697  if (fLineColor<0)
698  {
699   Float_t q=t->GetCharge();
700   curve->SetLineColor(kGreen);
701   if (q>0) curve->SetLineColor(kRed);
702   if (q<0) curve->SetLineColor(kBlue);
703  }
704  else
705  {
706   curve->SetLineColor(fLineColor);
707  }
708  curve->Draw();
709
710  if (!fCurves)
711  {
712   fCurves=new TObjArray();
713   fCurves->SetOwner();
714  }
715  fCurves->Add(curve);
716
717  // Display the marker for the track starting point
718  if (fMstyle>0)
719  {
720   TPolyMarker3D* m=new TPolyMarker3D();
721   m->SetPoint(0,points[0],points[1],points[2]);
722   m->SetMarkerStyle(fMstyle);
723   m->SetMarkerSize(fMsize);
724   Int_t col=curve->GetLineColor();
725   if (fMcol>0) col=fMcol;
726   m->SetMarkerColor(col);
727   m->Draw();
728   fCurves->Add(m);
729  }
730 }
731 ///////////////////////////////////////////////////////////////////////////
732 void AliHelix::Refresh(Int_t mode)
733 {
734 // Refresh the display screen before showing the next curve.
735 //
736 // mode = 0 : refreshing fully under user control.
737 //        1 : the display screen will be refreshed automatically
738 //            at each individual track display.
739 //       -1 : the display screen will be refreshed automatically
740 //            at each event display.
741 //
742 // The default is mode=0.
743
744  if (abs(mode)<2) fRefresh=mode;
745  if (fCurves) fCurves->Delete();
746 }
747 ///////////////////////////////////////////////////////////////////////////
748 void AliHelix::Display(AliEvent* evt,Double_t* range,Int_t iaxis,Double_t scale)
749 {
750 // Display the helix curves of all tracks of the specified event.
751 // Various events can be displayed together or individually; please refer to
752 // the memberfunction Refresh() for further details.
753 // Please refer to the track display memberfunction for further details
754 // on the input arguments.
755 // 
756 // The default values are range=0, iaxis=3 and scale=-1.
757 //
758 // Note :
759 // ------
760 // Before any display activity, a TCanvas and a TView have to be initiated
761 // first by the user like for instance
762 // 
763 // TCanvas* c1=new TCanvas("c1","c1");
764 // TView* view=new TView(1);
765 // view->SetRange(-1000,-1000,-1000,1000,1000,1000);
766 // view->ShowAxis();
767 //
768 // The user can also use the 3D viewing facilities from the TCanvas menu
769 // to open an appropriate view. 
770
771  if (!evt) return;
772
773  if (fRefresh<0) Refresh(fRefresh);
774
775  Int_t ntk=evt->GetNtracks();
776  for (Int_t jtk=1; jtk<=ntk; jtk++)
777  {
778   AliTrack* tx=evt->GetTrack(jtk);
779   if (tx) Display(tx,range,iaxis,scale);
780  }
781 }
782 ///////////////////////////////////////////////////////////////////////////
783 void AliHelix::Display(TObjArray* arr,Double_t* range,Int_t iaxis,Double_t scale)
784 {
785 // Display the helix curves of all tracks in the specified array.
786 // A convenient way to obtain an array with selected tracks from e.g. an AliEvent
787 // is to make use of its GetTracks() selection facility.
788 // Various arrays can be displayed together or individually; please refer to
789 // the memberfunction Refresh() for further details.
790 // Please refer to the track display memberfunction for further details
791 // on the input arguments.
792 // 
793 // The default values are range=0, iaxis=3 and scale=-1.
794 //
795 // Note :
796 // ------
797 // Before any display activity, a TCanvas and a TView have to be initiated
798 // first by the user like for instance
799 // 
800 // TCanvas* c1=new TCanvas("c1","c1");
801 // TView* view=new TView(1);
802 // view->SetRange(-1000,-1000,-1000,1000,1000,1000);
803 // view->ShowAxis();
804 //
805 // The user can also use the 3D viewing facilities from the TCanvas menu
806 // to open an appropriate view. 
807
808  if (!arr) return;
809
810  Int_t ntk=arr->GetEntries();
811  for (Int_t jtk=0; jtk<ntk; jtk++)
812  {
813   TObject* obj=arr->At(jtk);
814   if (!obj) continue;
815   if (!(obj->InheritsFrom("AliTrack"))) continue;
816   AliTrack* tx=(AliTrack*)obj;
817   Display(tx,range,iaxis,scale);
818  }
819 }
820 ///////////////////////////////////////////////////////////////////////////
821 AliPosition* AliHelix::Extrapolate(AliTrack* t,Double_t* pars,Double_t scale)
822 {
823 // Extrapolate an AliTrack according to the corresponding helix curve
824 // and provide a pointer to the impact position w.r.t. a specified plane.
825 // In case the track can never reach the specified plane, the returned
826 // position pointer is zero.
827 // Detailed information of all the helix points used in the extrapolation
828 // can be obtained via the GetN() and GetP() memberfunctions of TPolyLine3D.
829 // It is assumed that the track charge is stored in elementary units
830 // (i.e. charge=1 for a proton).
831 // The input argument "scale" specifies the unit scale for the various
832 // locations where scale=0.01 indicates unit scales in cm etc...
833 // In case scale<=0, the unit scale for locations is determined from the
834 // begin, reference or endpoint of the track. If neither of these
835 // positions is present, all locations are assumed to be given in meter.
836 // The extrapolation parameters for the impact plane and required accuracy
837 // are specified by pars[0], pars[1] and pars[2], respectively.
838 // pars[0] = coordinate value of the plane for the impact point
839 // pars[1] = required accuracy on the specified impact plane coordinate
840 // pars[2] = the axis along which the value of par[0] is specified
841 //
842 // The parameters can be specified either w.r.t. the LAB frame or the Helix frame.
843 // The latter is the frame in which the Z axis points in the B direction.
844 //
845 // The conventions for the par[2] argument are the following :
846 // par[2] = 1 ==> X axis in the LAB frame
847 //          2 ==> Y axis in the LAB frame
848 //          3 ==> Z axis in the LAB frame
849 //         -1 ==> X axis in the Helix frame
850 //         -2 ==> Y axis in the Helix frame
851 //         -3 ==> Z axis in the Helix frame
852 //
853 // Example :
854 // ---------
855 // To obtain an extrapolation to the plane Z=0 in the LAB frame
856 // with an accuracy of 0.001 cm the input arguments would be
857 // pars[0]=0  pars[1]=0.001  pars[2]=3  scale=0.01
858 //
859 // Note : The default value for the scale is -1.
860
861  if (fExt)
862  {
863   delete fExt;
864   fExt=0;
865  }
866
867  if (!t || !pars) return fExt;
868
869  AliPosition* rbeg=t->GetBeginPoint();
870  AliPosition* rend=t->GetEndPoint();
871  AliPosition* rref=t->GetReferencePoint();
872
873  // The unit scale for locations if not specified by the user
874  if (scale<=0)
875  {
876   scale=1; // Set default to meter
877   if (rbeg)
878   {
879    scale=rbeg->GetUnitScale();
880   }
881   else if (rend)
882   {
883    scale=rend->GetUnitScale();
884   }
885   else if (rref)
886   {
887    scale=rref->GetUnitScale();
888   }
889  }
890
891  Double_t range[2];
892  range[0]=pars[0]-fabs(pars[1])/2.;
893  range[1]=pars[0]+fabs(pars[1])/2.;
894
895  Int_t iaxis=int(pars[2]);
896
897  MakeCurve(t,range,iaxis,scale);
898
899  Int_t np=GetLastPoint()+1;
900  if (!np) return fExt;
901
902  Float_t* points=GetP();
903
904  // First point of the curve around the impact
905  Int_t ip=0;
906  Float_t first[3]={points[3*ip],points[3*ip+1],points[3*ip+2]};
907
908  // Last point of the curve around the impact
909  ip=GetLastPoint();
910  Float_t last[3]={points[3*ip],points[3*ip+1],points[3*ip+2]};
911
912  // The accuracy on the impact point 
913  Float_t err[3];
914  err[0]=fabs(first[0]-last[0]);
915  err[1]=fabs(first[1]-last[1]);
916  err[2]=fabs(first[2]-last[2]);
917
918  // Take the middle point as impact location
919  ip=np/2;
920  Float_t imp[3]={points[3*ip],points[3*ip+1],points[3*ip+2]};
921
922  fExt=new AliPosition();
923  fExt->SetUnitScale(scale);
924  fExt->SetPosition(imp,"car");
925  fExt->SetPositionErrors(err,"car");
926
927  return fExt;
928 }
929 ///////////////////////////////////////////////////////////////////////////