////
////
-
#include "AliTPCTracklet.h"
#include "TObjArray.h"
+#include "TLinearFitter.h"
#include "AliTPCseed.h"
#include "AliESDVertex.h"
+#include "AliTracker.h"
+#include "TTreeStream.h"
+#include "TRandom3.h"
+#include "TDecompChol.h"
+
+#include <iostream>
+using namespace std;
ClassImp(AliTPCTracklet)
+const Double_t AliTPCTracklet::kB2C=0.299792458e-3;
+
AliTPCTracklet::AliTPCTracklet()
- : fNClusters(0),fSector(-1),fOuter(0),fInner(0),fPrimary(0) {
+ : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(-1),fOuter(0),
+ fInner(0),fPrimary(0) {
////
// The default constructor. It is intended to be used for I/O only.
////
}
-AliTPCTracklet::AliTPCTracklet(const AliTPCseed *track,Int_t sector)
- : fNClusters(0),fSector(sector),fOuter(0),fInner(0),fPrimary(0) {
+AliTPCTracklet::AliTPCTracklet(const AliTPCseed *track,Int_t sector,
+ TrackType type,Bool_t storeClusters)
+ : fNClusters(0),fNStoredClusters(0),fClusters(0),fSector(sector),fOuter(0),
+ fInner(0),fPrimary(0) {
////
// Contructor for a tracklet out of a track. Only clusters within a given
// sector are used.
///
-
- AliTPCseed *t=new AliTPCseed(*track);
-
- if (!t->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-t->GetAlpha())) {
- delete t;
- return;
- }
- // fit from inner to outer row
- AliTPCseed *outerSeed=new AliTPCseed(*t);
- Int_t n=0;
+ //TODO: only kalman works
+
for (Int_t i=0;i<160;++i) {
- AliTPCclusterMI *c=t->GetClusterPointer(i);
- if (c&&c->GetDetector()==sector) {
- if (n==1) {
- outerSeed->ResetCovariance(100.);
- }
- ++n;
- Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
- Double_t cov[3]={0.1,0.,0.1}; //TODO: correct error parametrisation
- if (!outerSeed->PropagateTo(r[0])
- || !static_cast<AliExternalTrackParam*>(outerSeed)->Update(&r[1],cov)) {
- delete outerSeed;
- outerSeed=0;
- break;
- }
- }
+ AliTPCclusterMI *c=track->GetClusterPointer(i);
+ if (c&&c->GetDetector()==sector)
+ ++fNClusters;
}
- fNClusters=n;
- if (outerSeed)
- fOuter=new AliExternalTrackParam(*outerSeed);
- delete outerSeed;
- // fit from outer to inner rows
- AliTPCseed *innerSeed=new AliTPCseed(*t);
- n=0;
- for (Int_t i=159;i>=0;--i) {
- AliTPCclusterMI *c=t->GetClusterPointer(i);
- if (c&&c->GetDetector()==sector) {
- if (n==1) {
- innerSeed->ResetCovariance(100.);
- }
- ++n;
- Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
- Double_t cov[3]={0.1,0.,0.1};
- if (!innerSeed->PropagateTo(r[0])
- || !static_cast<AliExternalTrackParam*>(innerSeed)->Update(&r[1],cov)) {
- delete innerSeed;
- innerSeed=0;
- break;
- }
+
+ if (storeClusters) {
+ fClusters=new AliTPCclusterMI[fNClusters];
+ for (Int_t i=0;i<160;++i) {
+ AliTPCclusterMI *c=track->GetClusterPointer(i);
+ if (c&&c->GetDetector()==sector)
+ fClusters[fNStoredClusters]=c;
+ ++fNStoredClusters;
}
}
- fNClusters=TMath::Max(fNClusters,n);
- if (innerSeed)
- fInner=new AliExternalTrackParam(*outerSeed);
- // propagate to the primary vertex
- if (innerSeed) {
- AliTPCseed *primarySeed=new AliTPCseed(*innerSeed);
- Double_t pos[]={0.,0.,0.};
- Double_t sigma[]={.1,.1,.1}; //TODO: is this correct?
- AliESDVertex vertex(pos,sigma);
- if (primarySeed->PropagateToVertex(&vertex))
- fPrimary=new AliExternalTrackParam(*primarySeed);
- delete primarySeed;
+
+ switch (type) {
+ case kKalman:
+ FitKalman(track,sector);
+ break;
+ case kLinear:
+ case kQuadratic:
+ FitLinear(track,sector,type);
+ break;
+ case kRiemann:
+ FitRiemann(track,sector);
+ break;
}
- delete innerSeed;
- if (!fOuter&&!fInner)
- fNClusters=0;
+}
- delete t;
+AliTPCTracklet::AliTPCTracklet(const TObjArray &clusters,Int_t sector,
+ TrackType type,Bool_t storeClusters) {
+ //TODO: write it!
}
AliTPCTracklet::AliTPCTracklet(const AliTPCTracklet &t)
- : fNClusters(t.fNClusters),fSector(t.fSector),fOuter(0),fInner(0),
+ : fNClusters(t.fNClusters),fNStoredClusters(t.fNStoredClusters),fClusters(0),
+ fSector(t.fSector),fOuter(0),fInner(0),
fPrimary(0) {
////
// The copy constructor. You can copy tracklets!
////
+ if (t.fClusters) {
+ fClusters=new AliTPCclusterMI[t.fNStoredClusters];
+ for (int i=0;i<t.fNStoredClusters;++i)
+ fClusters[i]=t.fClusters[i];
+ }
if (t.fOuter)
fOuter=new AliExternalTrackParam(*t.fOuter);
if (t.fInner)
////
// The assignment constructor. You can assign tracklets!
////
- fNClusters=t.fNClusters;
- fSector=t.fSector;
if (this!=&t) {
+ fNClusters=t.fNClusters;
+ fNStoredClusters=fNStoredClusters;
+ delete fClusters;
+ if (t.fClusters) {
+ fClusters=new AliTPCclusterMI[t.fNStoredClusters];
+ for (int i=0;i<t.fNStoredClusters;++i)
+ fClusters[i]=t.fClusters[i];
+ }
+ else
+ fClusters=0;
+ fSector=t.fSector;
if (t.fOuter) {
if (fOuter)
*fOuter=*t.fOuter;
//
// The destructor. Yes, you can even destruct tracklets.
//
+ delete fClusters;
delete fOuter;
delete fInner;
delete fPrimary;
}
-TObjArray AliTPCTracklet::CreateTracklets(const AliTPCseed *s,
+void AliTPCTracklet::FitKalman(const AliTPCseed *track,Int_t sector) {
+ AliTPCseed *t=new AliTPCseed(*track);
+ if (!t->Rotate(TMath::DegToRad()*(sector%18*20.+10.)-t->GetAlpha())) {
+ delete t;
+ return;
+ }
+ // fit from inner to outer row
+ AliTPCseed *outerSeed=new AliTPCseed(*t);
+ Int_t n=0;
+ for (Int_t i=0;i<160;++i) {
+ AliTPCclusterMI *c=t->GetClusterPointer(i);
+ if (c&&c->GetDetector()==sector) {
+ if (n==1) {
+ outerSeed->ResetCovariance(100.);
+ }
+ ++n;
+ Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
+ Double_t cov[3]={0.1,0.,0.1}; //TODO: correct error parametrisation
+ if (!outerSeed->PropagateTo(r[0]) ||
+ !static_cast<AliExternalTrackParam*>(outerSeed)->Update(&r[1],cov)) {
+ delete outerSeed;
+ outerSeed=0;
+ break;
+ }
+ }
+ }
+ if (outerSeed)
+ fOuter=new AliExternalTrackParam(*outerSeed);
+ delete outerSeed;
+ // fit from outer to inner rows
+ AliTPCseed *innerSeed=new AliTPCseed(*t);
+ n=0;
+ for (Int_t i=159;i>=0;--i) {
+ AliTPCclusterMI *c=t->GetClusterPointer(i);
+ if (c&&c->GetDetector()==sector) {
+ if (n==1) {
+ innerSeed->ResetCovariance(100.);
+ }
+ ++n;
+ Double_t r[3]={c->GetX(),c->GetY(),c->GetZ()};
+ Double_t cov[3]={0.1,0.,0.1};
+ if (!innerSeed->PropagateTo(r[0]) ||
+ !static_cast<AliExternalTrackParam*>(innerSeed)->Update(&r[1],cov)) {
+ delete innerSeed;
+ innerSeed=0;
+ break;
+ }
+ }
+ }
+ if (innerSeed)
+ fInner=new AliExternalTrackParam(*innerSeed);
+ // propagate to the primary vertex
+ if (innerSeed) {
+ AliTPCseed *primarySeed=new AliTPCseed(*innerSeed);
+ Double_t pos[]={0.,0.,0.};
+ Double_t sigma[]={.1,.1,.1}; //TODO: is this correct?
+ AliESDVertex vertex(pos,sigma);
+ if (primarySeed->PropagateToVertex(&vertex))
+ fPrimary=new AliExternalTrackParam(*primarySeed);
+ delete primarySeed;
+ // for better comparison one does not want to have alpha changed...
+ if (!fPrimary->Rotate(fInner->GetAlpha())) {
+ delete fPrimary;
+ fPrimary=0;
+ }
+ }
+ delete innerSeed;
+
+ delete t;
+}
+
+void AliTPCTracklet::FitLinear(const AliTPCseed *track,Int_t sector,
+ TrackType type) {
+ TLinearFitter fy(1);
+ TLinearFitter fz(1);
+ fy.StoreData(kFALSE);
+ fz.StoreData(kFALSE);
+ switch (type) {
+ case kLinear:
+ fy.SetFormula("1 ++ x");
+ fz.SetFormula("1 ++ x");
+ break;
+ case kQuadratic:
+ fy.SetFormula("1 ++ x ++ x*x");
+ fz.SetFormula("1 ++ x");
+ break;
+ }
+ Double_t xmax=-1.;
+ Double_t xmin=1000.;
+ for (Int_t i=0;i<160;++i) {
+ AliTPCclusterMI *c=track->GetClusterPointer(i);
+ if (c&&c->GetDetector()==sector) {
+ Double_t x=c->GetX();
+ fy.AddPoint(&x,c->GetY());
+ fz.AddPoint(&x,c->GetZ());
+ xmax=TMath::Max(xmax,x);
+ xmin=TMath::Min(xmin,x);
+ }
+ }
+ fy.Eval();
+ fz.Eval();
+ Double_t a[3]={fy.GetParameter(0),
+ fy.GetParameter(1),
+ type==kQuadratic?fy.GetParameter(2):0.};
+ Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
+ fy.GetCovarianceMatrixElement(1,0),
+ fy.GetCovarianceMatrixElement(1,1),
+ type==kQuadratic?fy.GetCovarianceMatrixElement(2,0):0.,
+ type==kQuadratic?fy.GetCovarianceMatrixElement(2,1):0.,
+ type==kQuadratic?fy.GetCovarianceMatrixElement(2,2):0.};
+ for (int i=0;i<6;++i) ca[i]*=fy.GetChisquare()/fNClusters;
+ Double_t b[2]={fz.GetParameter(0),
+ fz.GetParameter(1)};
+ Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
+ fz.GetCovarianceMatrixElement(1,0),
+ fz.GetCovarianceMatrixElement(1,1)};
+ for (int i=0;i<3;++i) cb[i]*=fz.GetChisquare()/fNClusters;
+ Double_t p[5];
+ Double_t c[15];
+ Double_t alpha=track->GetAlpha();
+ Quadratic2Helix(a,ca,b,cb,0.,p,c);
+ fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
+ Quadratic2Helix(a,ca,b,cb,xmin,p,c);
+ fInner=new AliExternalTrackParam(xmin,alpha,p,c);
+ Quadratic2Helix(a,ca,b,cb,xmax,p,c);
+ fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
+}
+
+void AliTPCTracklet::Quadratic2Helix(Double_t *a,Double_t *ca,
+ Double_t *b,Double_t *cb,
+ Double_t x0,
+ Double_t *p,Double_t *c) {
+ // y(x)=a[0]+a[1]*x+a[2]*x^2
+ // z(x)=b[0]+b[1]*x
+ // parametrises the corosponding helix at x0
+
+ // get the polynoms at x0
+ Double_t a0=x0*x0*a[2] + x0*a[1] + a[0];
+ Double_t a1=2.*x0*a[2] + a[1];
+ Double_t a2= a[2];
+ Double_t ca00=ca[0]+x0*(2.*ca[1]+x0*(ca[2]+2.*ca[3]+x0*(2.*ca[4]+x0*ca[5])));
+ Double_t ca10=ca[1]+x0*(ca[2]+2.*ca[3]+x0*(3.*ca[4]+x0*2.*ca[5]));
+ Double_t ca11=ca[2]+x0*4.*(ca[4]+x0*ca[5]);
+ Double_t ca20=ca[3]+x0*(ca[4]+x0*ca[5]);
+ Double_t ca21=ca[3]+x0*2.*ca[5];
+ Double_t ca22=ca[5];
+
+ Double_t b0=x0*b[1] + b[0];
+ Double_t b1= b[1];
+ Double_t cb00=cb[0]+x0*(2.*cb[1]+x0*cb[2]);
+ Double_t cb10=cb[1]+x0*cb[2];
+ Double_t cb11=cb[2];
+
+ // transform to helix parameters
+ Double_t f =1.+a1*a1;
+ Double_t f2 =f*f;
+ Double_t fi =1./f;
+ Double_t fi12=TMath::Sqrt(fi);
+ Double_t fi32=fi*fi12;
+ Double_t fi2 =fi*fi;
+ Double_t fi52=fi2*fi12;
+ Double_t fi3 =fi2*fi;
+ Double_t fi5 =fi2*fi3;
+
+ Double_t xyz[3]={0.}; // TODO...
+ Double_t fc=1./(GetBz(xyz)*kB2C);
+
+ p[0]=a0; // y0
+ p[1]=b0; // z0
+ p[2]=a1*fi12; // snp
+ p[3]=b1; // tgl
+ p[4]=2.*a2*fi32*fc; // 1/pt
+
+ c[0] =ca00; // y0-y0
+ c[1] =0.; // z0-y0
+ c[2] =cb00; // z0-z0
+ c[3] =ca10*fi32; // snp-y0
+ c[4] =0.; // snp-z0
+ c[5] =ca11*fi3; // snp-snp
+ c[6] =0.; // tgl-y0
+ c[7] =cb10; // tgl-z0
+ c[8] =0.; // tgl-snp
+ c[9] =cb11; // tgl-tgl
+ c[10]=2.*(-3.*a1*a2*ca10+f*ca20)*fi3*fc; // 1/pt-y0
+ c[11]=0.; // 1/pt-z0
+ c[12]=2.*(-3.*a1*a2*ca11+f*ca21)*fi52*fc; // 1/pt-snp
+ c[13]=0.; // 1/pt-tgl
+ c[14]=(-12.*a1*a2*(-3.*a1*a2*ca11+2.*f*ca21)+4.*f2*ca22)*fi5
+ *fc*fc; // 1/pt-1/pt
+}
+
+
+void AliTPCTracklet::FitRiemann(const AliTPCseed *track,Int_t sector) {
+ TLinearFitter fy(2);
+ fy.StoreData(kFALSE);
+ fy.SetFormula("hyp2");
+ Double_t xmax=-1.;
+ Double_t xmin=1000.;
+ for (Int_t i=0;i<160;++i) {
+ AliTPCclusterMI *c=track->GetClusterPointer(i);
+ if (c&&c->GetDetector()==sector) {
+ Double_t x=c->GetX();
+ Double_t y=c->GetY();
+ Double_t xy[2]={x,y};
+ Double_t r=x*x+y*y;
+ Double_t errx=1.,erry=1.;//TODO!
+ Double_t err=TMath::Sqrt(4.*x*x*errx+4.*y*y*erry);
+ err=1.;
+ fy.AddPoint(xy,r,err);
+ xmax=TMath::Max(xmax,x);
+ xmin=TMath::Min(xmin,x);
+ }
+ }
+ fy.Eval();
+ Double_t a[3]={fy.GetParameter(0),
+ fy.GetParameter(1),
+ fy.GetParameter(2)};
+ Double_t ca[6]={fy.GetCovarianceMatrixElement(0,0),
+ fy.GetCovarianceMatrixElement(1,0),
+ fy.GetCovarianceMatrixElement(1,1),
+ fy.GetCovarianceMatrixElement(2,0),
+ fy.GetCovarianceMatrixElement(2,1),
+ fy.GetCovarianceMatrixElement(2,2)};
+
+ TLinearFitter fz(1);
+ fz.StoreData(kFALSE);
+ fz.SetFormula("hyp1");
+ Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
+ Double_t oldx=0.;
+ Double_t oldy=R;
+ Double_t phi=0.;
+ for (Int_t i=0;i<160;++i) {
+ AliTPCclusterMI *c=track->GetClusterPointer(i);
+ if (c&&c->GetDetector()==sector) {
+ Double_t x=c->GetX();
+ Double_t y=c->GetY();
+ Double_t dx=x-oldx;
+ Double_t dy=y-oldy;
+ phi+=2.*TMath::Abs(TMath::ATan2(.5*TMath::Sqrt(dx*dx+dy*dy),R));
+ Double_t err=1.;
+ fz.AddPoint(&phi,c->GetZ(),err);
+ oldx=x;
+ oldy=y;
+ }
+ }
+ fz.Eval();
+ Double_t b[2]={fz.GetParameter(0),
+ fz.GetParameter(1)};
+ Double_t cb[3]={fz.GetCovarianceMatrixElement(0,0),
+ fz.GetCovarianceMatrixElement(1,0),
+ fz.GetCovarianceMatrixElement(1,1)};
+
+ Double_t p[5];
+ Double_t c[15];
+ Double_t alpha=track->GetAlpha();
+ if (Riemann2Helix(a,ca,b,cb,0.,p,c))
+ fPrimary=new AliExternalTrackParam(0.,alpha,p,c);
+ if (Riemann2Helix(a,ca,b,cb,xmin,p,c))
+ fInner=new AliExternalTrackParam(xmin,alpha,p,c);
+ if (Riemann2Helix(a,ca,b,cb,xmax,p,c))
+ fOuter=new AliExternalTrackParam(xmax,alpha,p,c);
+}
+
+Bool_t AliTPCTracklet::Riemann2Helix(Double_t *a,Double_t *ca,
+ Double_t *b,Double_t *cb,
+ Double_t x0,
+ Double_t *p,Double_t *c) {
+ //TODO: signs!
+
+ Double_t xr0=.5*a[1];
+ Double_t yr0=.5*a[2];
+ Double_t R=.5*TMath::Sqrt(4.*a[0]+a[1]*a[1]+a[2]*a[2]);
+ Double_t dx=x0-xr0;
+ if (dx*dx>=R*R) return kFALSE;
+ Double_t dy=TMath::Sqrt(R*R-dx*dx); //sign!!
+ if (TMath::Abs(yr0+dy)>TMath::Abs(yr0-dy))
+ dy=-dy;
+ Double_t y0=yr0+dy;
+ Double_t tgp=-dx/dy; //TODO: dy!=0
+ Double_t z0=b[0]+TMath::ATan(tgp)*b[1];
+ Double_t xyz[3]={x0,y0,z0};
+ Double_t fc=1./(GetBz(xyz)*kB2C);
+ fc=1;
+ p[0]=y0; // y0
+ p[1]=z0; // z0
+ p[2]=tgp/TMath::Sqrt(1.+tgp*tgp); // snp
+ p[3]=b[1]; // tgl
+ p[4]=1./R*fc; // 1/pt
+
+ c[0] =0.; // y0-y0
+ c[1] =0.; // z0-y0
+ c[2] =0.; // z0-z0
+ c[3] =0.; // snp-y0
+ c[4] =0.; // snp-z0
+ c[5] =0.; // snp-snp
+ c[6] =0.; // tgl-y0
+ c[7] =0.; // tgl-z0
+ c[8] =0.; // tgl-snp
+ c[9] =0.; // tgl-tgl
+ c[10]=0.; // 1/pt-y0
+ c[11]=0.; // 1/pt-z0
+ c[12]=0.; // 1/pt-snp
+ c[13]=0.; // 1/pt-tgl
+ c[14]=0.; // 1/pt-1/pt
+
+ return kTRUE;
+}
+
+TObjArray AliTPCTracklet::CreateTracklets(const AliTPCseed *track,
+ TrackType type,
+ Bool_t storeClusters,
Int_t minClusters,
Int_t maxTracklets) {
// The tracklet factory: It creates several tracklets out of a track. They
Int_t sectors[72]={0};
for (Int_t i=0;i<160;++i) {
- AliTPCclusterMI *c=s->GetClusterPointer(i);
+ AliTPCclusterMI *c=track->GetClusterPointer(i);
if (c)
++sectors[c->GetDetector()];
}
TObjArray tracklets;
if (maxTracklets>72) maxTracklets=72; // just to protect against "users".
for (Int_t i=0;i<maxTracklets&§ors[indices[i]]>=minClusters;++i) {
- tracklets.Add(new AliTPCTracklet(s,indices[i]));
+ tracklets.Add(new AliTPCTracklet(track,indices[i],type,storeClusters));
}
return tracklets;
}
+
+TObjArray AliTPCTracklet::CreateTracklets(const TObjArray &clusters,
+ TrackType type,
+ Bool_t storeClusters,
+ Int_t minClusters,
+ Int_t maxTracklets) {
+
+}
+
+Bool_t AliTPCTracklet::PropagateToMeanX(const AliTPCTracklet &t1,
+ const AliTPCTracklet &t2,
+ AliExternalTrackParam *&t1m,
+ AliExternalTrackParam *&t2m) {
+ // This function propagates two Tracklets to a common x-coordinate. This
+ // x is dermined as the one that is in the middle of the two tracklets (they
+ // are assumed to live on two distinct x-intervalls).
+ // The inner parametrisation of the outer Tracklet and the outer
+ // parametrisation of the inner Tracklet are used and propagated to this
+ // common x. This result is saved not inside the Tracklets but two new
+ // ExternalTrackParams are created (that means you might want to delete
+ // them).
+ // In the case that the alpha angles of the Tracklets differ both angles
+ // are tried out for this propagation.
+ // In case of any failure kFALSE is returned, no AliExternalTrackParam
+ // is created und the pointers are set to 0.
+
+ if (t1.GetInner() && t1.GetOuter() &&
+ t2.GetInner() && t2.GetOuter()) {
+ if (t1.GetOuter()->GetX()<t2.GetInner()->GetX()) {
+ t1m=new AliExternalTrackParam(*t1.GetOuter());
+ t2m=new AliExternalTrackParam(*t2.GetInner());
+ }
+ else {
+ t1m=new AliExternalTrackParam(*t1.GetInner());
+ t2m=new AliExternalTrackParam(*t2.GetOuter());
+ }
+ Double_t mx=.5*(t1m->GetX()+t2m->GetX());
+ Double_t b1,b2;
+ Double_t xyz[3];
+ t1m->GetXYZ(xyz);
+ b1=GetBz(xyz);
+ t2m->GetXYZ(xyz);
+ b2=GetBz(xyz);
+ if (t1m->Rotate(t2m->GetAlpha())
+ && t1m->PropagateTo(mx,b1)
+ && t2m->PropagateTo(mx,b2));
+ else
+ if (t2m->Rotate(t1m->GetAlpha())
+ && t1m->PropagateTo(mx,b1)
+ && t2m->PropagateTo(mx,b2));
+ else {
+ delete t1m;
+ delete t2m;
+ t1m=t2m=0;
+ }
+ }
+ else {
+ t1m=t2m=0;
+ }
+ return t1m&&t2m;
+}
+
+double AliTPCTracklet::GetBz(Double_t *xyz) {
+ if (AliTracker::UniformField())
+ return AliTracker::GetBz();
+ else
+ return AliTracker::GetBz(xyz);
+}
+
+void AliTPCTracklet::RandomND(Int_t ndim,const Double_t *p,const Double_t *c,
+ Double_t *x) {
+ // This function generates a n-dimensional random variable x with mean
+ // p and covariance c.
+ // That is done using the cholesky decomposition of the covariance matrix,
+ // Begin_Latex C=U^{t} U End_Latex, with Begin_Latex U End_Latex being an
+ // upper triangular matrix. Given a vector v of iid gausian random variables
+ // with variance 1 one obtains the asked result as: Begin_Latex x=U^t v
+ // End_Latex.
+ // c is expected to be in a lower triangular format:
+ // c[0]
+ // c[1] c[2]
+ // c[3] c[4] c[5]
+ // etc.
+ static TRandom3 random;
+ Double_t *c2= new Double_t[ndim*ndim];
+ Int_t k=0;
+ for (Int_t i=0;i<ndim;++i)
+ for (Int_t j=0;j<=i;++j)
+ c2[i*ndim+j]=c2[j*ndim+i]=c[k++];
+ TMatrixDSym cm(ndim,c2);
+ delete[] c2;
+ TDecompChol chol(cm);
+ chol.Decompose();
+ const TVectorD pv(ndim);
+ const_cast<TVectorD*>(&pv)->Use(ndim,const_cast<Double_t*>(p));
+ TVectorD xv(ndim);
+ xv.Use(ndim,x);
+ for (Int_t i=0;i<ndim;++i)
+ xv[i]=random.Gaus();
+ TMatrixD L=chol.GetU();
+ L.T();
+ xv=L*xv+pv;
+}
+
+TEllipse AliTPCTracklet::ErrorEllipse(Double_t x,Double_t y,
+ Double_t sx,Double_t sy,Double_t sxy) {
+ /* Begin_Latex
+ r_{1,2}=1/2 (a+c#pm#sqrt{(a-c)^{2}+(2b)^{2}})
+ End_Latex */
+ Double_t det1=1./(sx*sy-sxy*sxy);
+ Double_t a=sy*det1;
+ Double_t b=-sxy*det1;
+ Double_t c=sx*det1;
+ Double_t d=c-a;
+ Double_t s=TMath::Sqrt(d*d+4.*b*b);
+ Double_t r1=TMath::Sqrt(.5*(a+c-s));
+ Double_t r2=TMath::Sqrt(.5*(a+c+s));
+ Double_t alpha=.5*TMath::ATan2(2.*b,d);
+ return TEllipse(x,y,r1,r2,0.,360.,alpha*TMath::RadToDeg());
+}
+
+void AliTPCTracklet::Test(const char* filename) {
+ /*
+ aliroot
+ AliTPCTracklet::Test("");
+ TFile f("AliTPCTrackletDebug.root");
+ TTree *t=f.Get("AliTPCTrackletDebug");
+ t->Draw("p0:p4");
+ TEllipse e=AliTPCTracklet::ErrorEllipse(0.,0.,4.,1.,1.8);
+ e.Draw();
+ */
+ TTreeSRedirector ds("AliTPCTrackletDebug.root");
+ Double_t p[5]={0.};
+ Double_t c[15]={4.,
+ 0.,4.,
+ 0.,0.,9.,
+ 0.,0.,0.,16.,
+ 1.8,0.,0.,0.,1.};
+ for (Int_t i=0;i<10000;++i) {
+ Double_t x[5];
+ RandomND(5,p,c,x);
+ ds<<"AliTPCTrackletDebug"
+ <<"p0="<<x[0]
+ <<"p1="<<x[1]
+ <<"p2="<<x[2]
+ <<"p3="<<x[3]
+ <<"p4="<<x[4]
+ <<"\n";
+ }
+
+ /*
+ Double_t b;
+ Double_t x=0.;
+ Double_t alpha=0.;
+ Double_t param[5]={0.};
+ Double_t covar[15]={1.,
+ 0.,1.,
+ 0.,0.,1.,
+ 0.,0.,0.,1.,
+ 0.,0.,0.,0.,1.};
+ AliExternalTrackParam track(x,alpha,param,covar);
+
+
+
+ for (Int_t i=0;i<points.GetNPoints();++i) {
+ Double_t x=0.;
+ Double_t alpha=0.;
+ Double_t param[5]={0.};
+ Double_t covar[15]={1.,
+ 0.,1.,
+ 0.,0.,1.,
+ 0.,0.,0.,1.,
+ 0.,0.,0.,0.,1.};
+ AliExternalTrackParam track(x,alpha,param,covar);
+ for (x=90.;x<250.;x+=1.) {
+ track.PropagateTo(x,b);
+ AliTPCclusterMI c();
+ }
+ }
+ */
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff