NcHelix.cxx

Go to the documentation of this file.

 
#include "NcHelix.h"
#include "Riostream.h"
 
ClassImp(NcHelix); // Class implementation to enable ROOT I/O
 
NcHelix::NcHelix() : THelix()
{
 
 fRefresh=0;
 fCurves=0;
 fExt=0;
 fTofmax=1e-8;
 fMstyle=-1;
 fMsize=0;
 fMcol=0;
 fEnduse=1;
 
 fLineWidth=2;
 fLineColor=-1;
}

NcHelix::~NcHelix()
{
 
 if (fCurves)
 {
  delete fCurves;
  fCurves=0;
 }
 if (fExt)
 {
  delete fExt;
  fExt=0;
 }
}

NcHelix::NcHelix(const NcHelix& h) : THelix(h)
{
 
 fB=h.fB;
 fRefresh=h.fRefresh;
 fTofmax=h.fTofmax;
 fMstyle=h.fMstyle;
 fMsize=h.fMsize;
 fMcol=h.fMcol;
 fEnduse=h.fEnduse;
}

void NcHelix::SetB(Nc3Vector& b)
{
 
 fB=b;
 
 if (fB.GetNorm()>0)
 {
  Double_t axis[3];
  fB.GetVector(axis,"car");
  SetAxis(axis);
 }
}

Nc3Vector& NcHelix::GetB()
{
 
 return fB;
}

void NcHelix::SetTofmax(Float_t tof)
{
 
 fTofmax=tof;
}

Float_t NcHelix::GetTofmax() const
{
 
 return fTofmax;
}

void NcHelix::SetMarker(Int_t style,Float_t size,Int_t col)
{
 
 fMstyle=style;
 fMsize=size;
 fMcol=col;
}

void NcHelix::UseEndPoint(Int_t mode)
{
 
 if (mode==0 || mode==1) fEnduse=mode; 
}

void NcHelix::MakeCurve(NcTrack* t,Double_t* range,Int_t iaxis,Double_t scale)
{
 
 SetPolyLine(0); // Reset the polyline data points
 
 if (!t || (range && !iaxis)) return;
 
 Double_t energy=t->GetEnergy(1); // Track energy in GeV
 Double_t betanorm=t->GetBeta();
 
 if (energy<=0 || betanorm<=0) return;
 
 NcPosition* rbeg=t->GetBeginPoint();
 NcPosition* rend=0;
 if (fEnduse) rend=t->GetEndPoint();
 NcPosition* rref=t->GetReferencePoint();
 
 // Magnetic field vector or default Z-direction
 Double_t bvec[3]={0,0,1};
 if (fB.GetNorm()>0) fB.GetVector(bvec,"car");
 
 // The unit scale for locations if not specified by the user
 if (scale<=0)
 {
  scale=1; // Set default to meter
  if (rbeg)
  {
   scale=rbeg->GetUnitScale();
  }
  else if (rend)
  {
   scale=rend->GetUnitScale();
  }
  else if (rref)
  {
   scale=rref->GetUnitScale();
  }
 }
 
 Double_t c=2.99792458e8/scale; // Lightspeed in the selected unit scale
 
 // The helix angular frequency
 Double_t w=9e7*(t->GetCharge()*fB.GetNorm())/energy;
 
 // The particle velocity in the LAB frame
 Nc3Vector beta=t->GetBetaVector();
 Nc3Vector v=beta*c;
 Double_t vel[3];
 v.GetVector(vel,"car");
 
 // The particle velocity in the Helix frame
 Nc3Vector betaprim=beta.GetPrimed(fRotMat);
 v=v.GetPrimed(fRotMat);
 Double_t velprim[3];
 v.GetVector(velprim,"car");
 
 // Check compatibility of velocity and range specification.
 if (range)
 {
  Double_t betavec[3];
  if (iaxis>0) beta.GetVector(betavec,"car");
  if (iaxis<0) betaprim.GetVector(betavec,"car");
  if (fabs(betavec[abs(iaxis)-1])/betanorm<1e-10) return;
 }
 
 // The LAB location in which the velocity of the particle is defined
 Double_t loc[3]={0,0,0};
 Nc3Vector rx;
 Double_t scalex=-1;
 if (rref)
 {
  rx=(Nc3Vector)(*rref);
  scalex=rref->GetUnitScale();
 }
 else if (rbeg)
 {
  rx=(Nc3Vector)(*rbeg);
  scalex=rbeg->GetUnitScale();
 }
 else if (rend)
 {
  rx=(Nc3Vector)(*rend);
  scalex=rend->GetUnitScale();
 }
 
 if (scalex>0 && (scalex/scale>1.1 || scale/scalex>1.1)) rx*=scalex/scale;
 rx.GetVector(loc,"car");
 
 // Initialisation of Helix kinematics
 SetHelix(loc,vel,w,0,kUnchanged,bvec);
 
 Int_t bend=0;
 if (fabs(w)>0 && fabs(fVt)>0) bend=1;
 
 // Flight time boundaries.
 // The time origin t=0 is chosen to indicate the position in which
 // the particle velocity was defined.
 // The total flight time is initialised to the (user specified) tofmax.
 Double_t tmin=0,tmax=0;
 Double_t tof=fTofmax;
 Double_t dum=0;
 
 // The trajectory begin and end points
 Double_t vec1[3]={0,0,0};
 Double_t vec2[3]={0,0,0};
 Nc3Vector r1;
 Nc3Vector r2;
 Double_t scale1=1;
 Double_t scale2=1;
 
 if (!bend)
 {
  // Treatment of straight trajectories //
  Nc3Vector r;
  if (range) // Specified range allows for exact flight time boundaries
  {
   if (iaxis>0)
   {
    tmin=(range[0]-loc[iaxis-1])/vel[iaxis-1];
    tmax=(range[1]-loc[iaxis-1])/vel[iaxis-1];
   }
   else
   {
    loc[0]=fX0;
    loc[1]=fY0;
    loc[2]=fZ0;
    tmin=(range[0]-loc[abs(iaxis)-1])/velprim[abs(iaxis)-1];
    tmax=(range[1]-loc[abs(iaxis)-1])/velprim[abs(iaxis)-1];
   }
   if (tmax<tmin)
   {
    dum=tmin;
    tmin=tmax;
    tmax=dum;
   }
   // Make the 'curve' in the LAB frame and exit.
   // Use the parametrisation : r(t)=r0+t*v
   // using the range based flight time boundaries.
   // An additional point in the middle of the trajectory is
   // generated in view of accuracy in the case of extrapolations.
   tof=tmax-tmin;
   v=beta*c;
   r1=rx;
   r=v*tmin;
   r1=r1+r;
   r1.GetVector(vec1,"car");
   SetNextPoint(float(vec1[0]),float(vec1[1]),float(vec1[2]));
   r=v*(tof/2.);
   r2=r1+r;
   r2.GetVector(vec2,"car");
   SetNextPoint(float(vec2[0]),float(vec2[1]),float(vec2[2]));
   r=v*tof;
   r2=r1+r;
   r2.GetVector(vec2,"car");
   SetNextPoint(float(vec2[0]),float(vec2[1]),float(vec2[2]));
  }
  else // Automatic range determination
  {
   // Initially the point with Z=0 in the Helix frame is taken as a starting point.
   // In case this point can't be reached, the point in which the particle velocity
   // was defined is taken as the starting point.
   // The endpoint is initially obtained by applying the tofmax from the start point.
   tmin=0;
   if (fabs(fVz)>0) tmin=-fZ0/fVz;
   v=beta*c;
   r1=rx;
   r=v*tmin;
   r1=r1+r;
 
   // Override the initial begin and endpoint settings by the track data
   if (rbeg)
   {
    r1=(Nc3Vector)(*rbeg); 
    scale1=rbeg->GetUnitScale();
    // All coordinates in the selected unit scale
    if (scale1/scale>1.1 || scale/scale1>1.1) r1*=scale1/scale;
   }
 
   r=v*fTofmax;
   r2=r1+r;
   if (rend)
   {
    r2=(Nc3Vector)(*rend); 
    scale2=rend->GetUnitScale();
    // All coordinates in the selected unit scale
    if (scale2/scale>1.1 || scale/scale2>1.1) r2*=scale2/scale;
   }
   
   r1.GetVector(vec1,"car");
   r2.GetVector(vec2,"car");
 
   // Make the 'curve' in the LAB frame and exit.
   SetNextPoint(float(vec1[0]),float(vec1[1]),float(vec1[2]));
   SetNextPoint(float(vec2[0]),float(vec2[1]),float(vec2[2]));
  }
 }
 else
 {
  // Treatment of curved trajectories //
 
  // Initialisation of the flight time boundaries.
  // Based on the constant motion of the particle along the Helix Z-axis,
  // the parametrisation z(t)=z0+fVz*t in the Helix frame is used. 
  // If possible the point with Z=0 in the Helix frame is taken as a starting point.
  // In case this point can't be reached, the point in which the particle velocity
  // was defined is taken as the starting point.
  tmin=0;
  if (fabs(fVz)>0) tmin=-fZ0/fVz;
  tmax=tmin+fTofmax;
 
  if (tmax<tmin)
  {
   dum=tmin;
   tmin=tmax;
   tmax=dum;
  }
 
  // Determination of the range in the helix frame
 
  if (!range) // Automatic range determination
  {
   scale1=1;
   scale2=1;
   if (rbeg)
   {
    r1=rbeg->GetPrimed(fRotMat);
    scale1=rbeg->GetUnitScale();
    // All coordinates in the selected unit scale
    if (scale1/scale>1.1 || scale/scale1>1.1) r1*=scale1/scale;
    // Re-calculate the tmin for this new starting point
    r1.GetVector(vec1,"car");
    if (fabs(fVz)>0) tmin=(vec1[2]-fZ0)/fVz;
    tmax=tmin+fTofmax;
   }
   if (rend)
   {
    r2=rend->GetPrimed(fRotMat);
    scale2=rend->GetUnitScale();
    // All coordinates in the selected unit scale
    if (scale2/scale>1.1 || scale/scale2>1.1) r2*=scale2/scale;
    r2.GetVector(vec2,"car");
    if (fabs(fVz)>0) tmax=(vec2[2]-fZ0)/fVz;
   }
   // Make the curve on basis of the flight time boundaries and exit
   if (tmax<tmin)
   {
    dum=tmin;
    tmin=tmax;
    tmax=dum;
   }
   SetRange(tmin,tmax,kHelixT);
  }
  else // User explicitly specified range
  {
   vec1[abs(iaxis)-1]=range[0];
   vec2[abs(iaxis)-1]=range[1];
   r1.SetVector(vec1,"car");
   r2.SetVector(vec2,"car");
   if (iaxis>0) // Range specified in LAB frame
   {
    r1=r1.GetPrimed(fRotMat);
    r1.GetVector(vec1,"car");
    r2=r2.GetPrimed(fRotMat);
    r2.GetVector(vec2,"car");
   } 
   // Determination of the axis component with the
   // largest range difference
   Double_t dmax=0;
   Int_t imax=0;
   Double_t test=0;
   for (Int_t i=0; i<3; i++)
   {
    test=fabs(vec1[i]-vec2[i]);
    if (test>dmax)
    {
     dmax=test;
     imax=i;
    }
   }
 
   Double_t rmin=vec1[imax];
   Double_t rmax=vec2[imax];
   if (rmax<rmin)
   {
    dum=rmin;
    rmin=rmax;
    rmax=dum;
   }
 
   // The kinematic range boundaries in the helix frame
   Double_t xmin=fX0-fVt/fW;
   Double_t xmax=fX0+fVt/fW;
   Double_t ymin=fY0-fVt/fW;
   Double_t ymax=fY0+fVt/fW;
 
   if (xmax<xmin)
   {
    dum=xmin;
    xmin=xmax;
    xmax=dum;
   }
   if (ymax<ymin)
   {
    dum=ymin;
    ymin=ymax;
    ymax=dum;
   }
 
   // Set the range for the helix
   if (imax==2 && dmax>0) SetRange(rmin,rmax,kHelixZ);
   if (imax==1)
   {
    // Limit range to kinematic boundaries if needed
    if (rmin<=ymin) rmin=ymin+1e-6*dmax;
    if (rmax>=ymax) rmax=ymax-1e-6*dmax;
    if (rmin<rmax) SetRange(rmin,rmax,kHelixY);
   }
   if (imax==0)
   {
    // Limit range to kinematic boundaries if needed
    if (rmin<=xmin) rmin=xmin+1e-6*dmax;
    if (rmax>=xmax) rmax=xmax-1e-6*dmax;
    if (rmin<rmax) SetRange(rmin,rmax,kHelixX);
   }
  }
 }
 return;
}

void NcHelix::Display(NcTrack* t,Double_t* range,Int_t iaxis,Double_t scale)
{
 
 if (!t || (range && !iaxis)) return;
 
 MakeCurve(t,range,iaxis,scale);
 
 if (fRefresh>0) Refresh(fRefresh);
 
 Int_t np=GetLastPoint()+1;
 if (!np) return;
 
 Float_t* points=GetP();
 TPolyLine3D* curve=new TPolyLine3D(np,points);
 
 curve->SetLineWidth(fLineWidth);
 if (fLineColor<0)
 {
  Float_t q=t->GetCharge();
  curve->SetLineColor(kGreen);
  if (q>0) curve->SetLineColor(kRed);
  if (q<0) curve->SetLineColor(kBlue);
 }
 else
 {
  curve->SetLineColor(fLineColor);
 }
 curve->Draw();
 
 if (!fCurves)
 {
  fCurves=new TObjArray();
  fCurves->SetOwner();
 }
 fCurves->Add(curve);
 
 // Display the marker for the track starting point
 if (fMstyle>0)
 {
  TPolyMarker3D* m=new TPolyMarker3D();
  m->SetPoint(0,points[0],points[1],points[2]);
  m->SetMarkerStyle(fMstyle);
  m->SetMarkerSize(fMsize);
  Int_t col=curve->GetLineColor();
  if (fMcol>0) col=fMcol;
  m->SetMarkerColor(col);
  m->Draw();
  fCurves->Add(m);
 }
}

void NcHelix::Refresh(Int_t mode)
{
 
 if (abs(mode)<2) fRefresh=mode;
 if (fCurves) fCurves->Delete();
}

void NcHelix::Display(NcEvent* evt,Double_t* range,Int_t iaxis,Double_t scale)
{
 
 if (!evt) return;
 
 if (fRefresh<0) Refresh(fRefresh);
 
 Int_t ntk=evt->GetNtracks();
 for (Int_t jtk=1; jtk<=ntk; jtk++)
 {
  NcTrack* tx=evt->GetTrack(jtk);
  if (tx) Display(tx,range,iaxis,scale);
 }
}

void NcHelix::Display(TObjArray* arr,Double_t* range,Int_t iaxis,Double_t scale)
{
 
 if (!arr) return;
 
 Int_t ntk=arr->GetEntries();
 for (Int_t jtk=0; jtk<ntk; jtk++)
 {
  TObject* obj=arr->At(jtk);
  if (!obj) continue;
  if (!(obj->InheritsFrom("NcTrack"))) continue;
  NcTrack* tx=(NcTrack*)obj;
  Display(tx,range,iaxis,scale);
 }
}

NcPosition* NcHelix::Extrapolate(NcTrack* t,Double_t* pars,Double_t scale)
{
 
 if (fExt)
 {
  delete fExt;
  fExt=0;
 }
 
 if (!t || !pars) return fExt;
 
 NcPosition* rbeg=t->GetBeginPoint();
 NcPosition* rend=t->GetEndPoint();
 NcPosition* rref=t->GetReferencePoint();
 
 // The unit scale for locations if not specified by the user
 if (scale<=0)
 {
  scale=1; // Set default to meter
  if (rbeg)
  {
   scale=rbeg->GetUnitScale();
  }
  else if (rend)
  {
   scale=rend->GetUnitScale();
  }
  else if (rref)
  {
   scale=rref->GetUnitScale();
  }
 }
 
 Double_t range[2];
 range[0]=pars[0]-fabs(pars[1])/2.;
 range[1]=pars[0]+fabs(pars[1])/2.;
 
 Int_t iaxis=int(pars[2]);
 
 MakeCurve(t,range,iaxis,scale);
 
 Int_t np=GetLastPoint()+1;
 if (!np) return fExt;
 
 Float_t* points=GetP();
 
 // First point of the curve around the impact
 Int_t ip=0;
 Float_t first[3]={points[3*ip],points[3*ip+1],points[3*ip+2]};
 
 // Last point of the curve around the impact
 ip=GetLastPoint();
 Float_t last[3]={points[3*ip],points[3*ip+1],points[3*ip+2]};
 
 // The accuracy on the impact point 
 Float_t err[3];
 err[0]=fabs(first[0]-last[0]);
 err[1]=fabs(first[1]-last[1]);
 err[2]=fabs(first[2]-last[2]);
 
 // Take the middle point as impact location
 ip=np/2;
 Float_t imp[3]={points[3*ip],points[3*ip+1],points[3*ip+2]};
 
 fExt=new NcPosition();
 fExt->SetUnitScale(scale);
 fExt->SetPosition(imp,"car");
 fExt->SetPositionErrors(err,"car");
 
 return fExt;
}