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