[u/mrichter/AliRoot.git] / RICH / AliRICHHelix.h

#ifndef AliRICHHelix_h
#define AliRICHHelix_h

#include <TObject.h>              //base class
#include <TVector3.h>             //used extensively
#include "AliRICHParam.h"         //RichIntersect() 

class AliRICHHelix: public TObject
{
public:
  AliRICHHelix():TObject()                                                                                                                  {}
  AliRICHHelix(Double_t p,Double_t theta,Double_t phi,Double_t bz=0.2);//p [GeV], theta,phi [deg], Bz [Tesla];
  AliRICHHelix(const TVector3 &x0,const TVector3 &p0,Int_t q=1,Double_t b=0.2):TObject(),fX0(x0),fP0(p0),fX(x0),fP(p0),fLen(0),fQ(q),fBz(b) {}
  virtual ~AliRICHHelix()                                                                                                                   {}        
           
         void     Draw          (const Option_t *opt=""              );                              //from TObject, draw helix
         void     Print         (const Option_t *opt=""              )const;                         //from TObject, print status
//private part++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  inline Bool_t   Intersection  (TVector3 planePnt,TVector3 planeNorm);                              //intersection with plane given by point and normal vector
  inline void     Propagate     (Double_t len                        );                              //propogate helix by given length along it
  inline Int_t    RichIntersect (AliRICHParam *pParam                );                              //search intersection with any RICH chamber 
         TVector2 PosRad        (                                    )const{return fPosRad;}         //intersection with radiator LORS
         TVector2 PosPc         (                                    )const{return fPosPc;}          //returns position of intersection with PC (local system)
         TVector3 Ploc          (                                    )const{return fPloc;}           //returns momentum at the position of intersection with radiator   
//         Double_t Length        ()                const{return fLen;}             //returns length of the track from initial point to RICH
         TVector3 X             (                                    )const{return fX;}
protected:
  TVector3 fX0;            //helix position in point of definition, [cm]    in MARS
  TVector3 fP0;            //helix momentum in point of definition, [GeV/c] in MARS
  TVector3 fX;             //helix position in point of interest,  [cm]    in MARS
  TVector3 fP;             //helix momentum in point of interest,  [GeV/c] in MARS
  Double_t fLen;           //helix length from point of definition to point of interest, [cm]    
  Int_t    fQ;             //sign of track charge (value not provided by current ESD)  
  Double_t fBz;            //magnetic field along z, [Tesla] 
  TVector2 fPosRad;        //helix intersection with radiator entrance, LORS [cm]
  TVector2 fPosPc;         //helix intersection with PC, LORS [cm]
  TVector3 fPloc;          //helix momentum, LORS [GeV/c]
  ClassDef(AliRICHHelix,0) //General helix
};//class AliRICHHelix
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
void AliRICHHelix::Propagate(Double_t len)
{
// Propogates the helix from inintial point by a given distance along helix. Assumes uniform magnetic field along z direction.  
// Arguments: len - distance to propagate by, [cm]
//   Returns: none    
  if(fBz==0){//no magnetic field->straight line
    fX=fX0+fP0.Unit()*len;   
  }else{
    const Double_t c = 0.00299792458;//this speed of light value provides that coordinates are in cm momentum in GeV/c
    Double_t a = -c*fBz*fQ;
     Double_t rho = a/fP0.Mag();
    fX.SetX( fX0.X()+fP0.X()*TMath::Sin(rho*len)/a-fP0.Y()*(1-TMath::Cos(rho*len))/a  );
    fX.SetY( fX0.Y()+fP0.Y()*TMath::Sin(rho*len)/a+fP0.X()*(1-TMath::Cos(rho*len))/a  ); 
    fX.SetZ( fX0.Z()+fP0.Z()*len/fP0.Mag()                                            );
    fP.SetX( fP0.X()*TMath::Cos(rho*len)-fP0.Y()*TMath::Sin(rho*len)                  );
    fP.SetY( fP0.Y()*TMath::Cos(rho*len)+fP0.X()*TMath::Sin(rho*len)                  );
    fP.SetZ( fP0.Z()                                                                  );
    fLen=len;
  }
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Bool_t AliRICHHelix::Intersection(TVector3 planePoint,TVector3 planeNorm)
{
// Finds point of intersection (if exists) of the helix with the plane. Stores result in fX and fP.   
// Arguments: planePoint,planeNorm - the plane defined by any plane's point and vector, normal to the plane
//   Returns: kTrue if helix intersects the plane, kFALSE otherwise.
  
  Double_t s=(planePoint-fX0)*planeNorm,dist=99999,distPrev=dist;//estimates initial distance to plane

  while(TMath::Abs(dist)>0.00001){
    Propagate(s);                        //calculates helix at the distance s from x0 ALONG the helix
    dist=(fX-planePoint)*planeNorm;      //distance between current helix position and plane
    if(TMath::Abs(dist) >= TMath::Abs(distPrev)) { return kFALSE;}
    distPrev=dist;
    s-=dist;
  }
  return kTRUE;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Int_t AliRICHHelix::RichIntersect(AliRICHParam *pParam)
{
// Searchs for intersection of this helix with all RICH chambers, returns chamber number or 0 if no intersection
// On exit fPosRad contain position of intersection in radiator LORS (cm)    
//         fPosPc  contains the same for photocathode  
  for(Int_t iCh=1;iCh<=AliRICHParam::kNch;iCh++){//chambers loop
    TVector3 norm =pParam->Lors2MarsVec(iCh,TVector3(0,0,1));
    if(Intersection(pParam->Lors2Mars(iCh,0,0,AliRICHParam::kRad),norm)){//there is intersection with radiator plane
      fPosRad=pParam->Mars2Lors(iCh,fX,AliRICHParam::kRad);//position on radiator plane
      if(pParam->IsAccepted(fPosRad)){//intersection within radiator (even if in dead zone)
        
        if(Intersection(pParam->Lors2Mars(iCh,0,0,AliRICHParam::kPc),norm)){//there is intersection with photocathode
          fPosPc=pParam->Mars2Lors(iCh,fX,AliRICHParam::kPc);//position on radiator plane
          if(pParam->IsAccepted(fPosPc)){//intersection within pc (even if in dead zone)
            
            fPloc=pParam->Mars2LorsVec(iCh,fP);//trasform p to local system
            return iCh;
          }//if inside PC
        }//if intersects PC
        
      }//if inside radiator
    }//if for radiator       
  }//chambers loop
  return 0;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#endif
Commit	Line	Data
998b831f	1	#ifndef AliRICHHelix_h
	2	#define AliRICHHelix_h
	3
db910db9	4	#include <TObject.h> //base class
	5	#include <TVector3.h> //used extensively
	6	#include "AliRICHParam.h" //RichIntersect()
998b831f	7
	8	class AliRICHHelix: public TObject
	9	{
	10	public:
db910db9	11	AliRICHHelix():TObject() {}
	12	AliRICHHelix(Double_t p,Double_t theta,Double_t phi,Double_t bz=0.2);//p [GeV], theta,phi [deg], Bz [Tesla];
	13	AliRICHHelix(const TVector3 &x0,const TVector3 &p0,Int_t q=1,Double_t b=0.2):TObject(),fX0(x0),fP0(p0),fX(x0),fP(p0),fLen(0),fQ(q),fBz(b) {}
	14	virtual ~AliRICHHelix() {}
d3eb6079	15
db910db9	16	void Draw (const Option_t *opt="" ); //from TObject, draw helix
	17	void Print (const Option_t *opt="" )const; //from TObject, print status
	18	//private part++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	19	inline Bool_t Intersection (TVector3 planePnt,TVector3 planeNorm); //intersection with plane given by point and normal vector
	20	inline void Propagate (Double_t len ); //propogate helix by given length along it
	21	inline Int_t RichIntersect (AliRICHParam *pParam ); //search intersection with any RICH chamber
	22	TVector2 PosRad ( )const{return fPosRad;} //intersection with radiator LORS
	23	TVector2 PosPc ( )const{return fPosPc;} //returns position of intersection with PC (local system)
	24	TVector3 Ploc ( )const{return fPloc;} //returns momentum at the position of intersection with radiator
	25	// Double_t Length () const{return fLen;} //returns length of the track from initial point to RICH
	26	TVector3 X ( )const{return fX;}
998b831f	27	protected:
db910db9	28	TVector3 fX0; //helix position in point of definition, [cm] in MARS
	29	TVector3 fP0; //helix momentum in point of definition, [GeV/c] in MARS
	30	TVector3 fX; //helix position in point of interest, [cm] in MARS
	31	TVector3 fP; //helix momentum in point of interest, [GeV/c] in MARS
	32	Double_t fLen; //helix length from point of definition to point of interest, [cm]
998b831f	33	Int_t fQ; //sign of track charge (value not provided by current ESD)
db910db9	34	Double_t fBz; //magnetic field along z, [Tesla]
	35	TVector2 fPosRad; //helix intersection with radiator entrance, LORS [cm]
	36	TVector2 fPosPc; //helix intersection with PC, LORS [cm]
	37	TVector3 fPloc; //helix momentum, LORS [GeV/c]
998b831f	38	ClassDef(AliRICHHelix,0) //General helix
998b831f	39	};//class AliRICHHelix
db910db9	40	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
d3eb6079	41	void AliRICHHelix::Propagate(Double_t len)
998b831f	42	{
db910db9	43	// Propogates the helix from inintial point by a given distance along helix. Assumes uniform magnetic field along z direction.
	44	// Arguments: len - distance to propagate by, [cm]
	45	// Returns: none
	46	if(fBz==0){//no magnetic field->straight line
	47	fX=fX0+fP0.Unit()*len;
	48	}else{
	49	const Double_t c = 0.00299792458;//this speed of light value provides that coordinates are in cm momentum in GeV/c
	50	Double_t a = -cfBzfQ;
	51	Double_t rho = a/fP0.Mag();
	52	fX.SetX( fX0.X()+fP0.X()TMath::Sin(rholen)/a-fP0.Y()(1-TMath::Cos(rholen))/a );
	53	fX.SetY( fX0.Y()+fP0.Y()TMath::Sin(rholen)/a+fP0.X()(1-TMath::Cos(rholen))/a );
	54	fX.SetZ( fX0.Z()+fP0.Z()*len/fP0.Mag() );
	55	fP.SetX( fP0.X()TMath::Cos(rholen)-fP0.Y()TMath::Sin(rholen) );
	56	fP.SetY( fP0.Y()TMath::Cos(rholen)+fP0.X()TMath::Sin(rholen) );
	57	fP.SetZ( fP0.Z() );
	58	fLen=len;
	59	}
998b831f	60	}
db910db9	61	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
998b831f	62	Bool_t AliRICHHelix::Intersection(TVector3 planePoint,TVector3 planeNorm)
998b831f	63	{
db910db9	64	// Finds point of intersection (if exists) of the helix with the plane. Stores result in fX and fP.
	65	// Arguments: planePoint,planeNorm - the plane defined by any plane's point and vector, normal to the plane
	66	// Returns: kTrue if helix intersects the plane, kFALSE otherwise.
998b831f	67
d3eb6079	68	Double_t s=(planePoint-fX0)*planeNorm,dist=99999,distPrev=dist;//estimates initial distance to plane
998b831f	69
d3eb6079	70	while(TMath::Abs(dist)>0.00001){
998b831f	71	Propagate(s); //calculates helix at the distance s from x0 ALONG the helix
998b831f	72	dist=(fX-planePoint)*planeNorm; //distance between current helix position and plane
e84e7c82	73	if(TMath::Abs(dist) >= TMath::Abs(distPrev)) { return kFALSE;}
998b831f	74	distPrev=dist;
	75	s-=dist;
	76	}
	77	return kTRUE;
	78	}
db910db9	79	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
998b831f	80	Int_t AliRICHHelix::RichIntersect(AliRICHParam *pParam)
	81	{
	82	// Searchs for intersection of this helix with all RICH chambers, returns chamber number or 0 if no intersection
db910db9	83	// On exit fPosRad contain position of intersection in radiator LORS (cm)
998b831f	84	// fPosPc contains the same for photocathode
cb801d41	85	for(Int_t iCh=1;iCh<=AliRICHParam::kNch;iCh++){//chambers loop
	86	TVector3 norm =pParam->Lors2MarsVec(iCh,TVector3(0,0,1));
	87	if(Intersection(pParam->Lors2Mars(iCh,0,0,AliRICHParam::kRad),norm)){//there is intersection with radiator plane
	88	fPosRad=pParam->Mars2Lors(iCh,fX,AliRICHParam::kRad);//position on radiator plane
998b831f	89	if(pParam->IsAccepted(fPosRad)){//intersection within radiator (even if in dead zone)
d3eb6079	90
cb801d41	91	if(Intersection(pParam->Lors2Mars(iCh,0,0,AliRICHParam::kPc),norm)){//there is intersection with photocathode
cb801d41	92	fPosPc=pParam->Mars2Lors(iCh,fX,AliRICHParam::kPc);//position on radiator plane
998b831f	93	if(pParam->IsAccepted(fPosPc)){//intersection within pc (even if in dead zone)
d3eb6079	94
cb801d41	95	fPloc=pParam->Mars2LorsVec(iCh,fP);//trasform p to local system
cb801d41	96	return iCh;
998b831f	97	}//if inside PC
d3eb6079	98	}//if intersects PC
d3eb6079	99
998b831f	100	}//if inside radiator
998b831f	101	}//if for radiator
cb801d41	102	}//chambers loop
998b831f	103	return 0;
998b831f	104	}
db910db9	105	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
998b831f	106	#endif