Adding new classes for TPC alingment - Magnus Mager
authormarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 1 Feb 2008 18:29:56 +0000 (18:29 +0000)
committermarian <marian@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 1 Feb 2008 18:29:56 +0000 (18:29 +0000)
TPC/AliTPCTracklet.cxx
TPC/AliTPCTracklet.h
TPC/AliTPCcalibAlign.cxx [new file with mode: 0644]
TPC/AliTPCcalibAlign.h [new file with mode: 0644]
TPC/AliTPCcalibAlignment.cxx [new file with mode: 0644]
TPC/AliTPCcalibAlignment.h [new file with mode: 0644]
TPC/TPCcalibLinkDef.h
TPC/libTPCcalib.pkg

index 0cef3d36e4d905287efd96ee03213e4a12fae945..7b3a523b7a1036990e4a5a2d894bc1ded978180d 100755 (executable)
 ////
 //// 
 
-
 #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)
@@ -134,9 +118,18 @@ AliTPCTracklet& AliTPCTracklet::operator=(const AliTPCTracklet &t) {
   ////
   // 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;
@@ -177,12 +170,323 @@ AliTPCTracklet::~AliTPCTracklet() {
   //
   // 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
@@ -194,7 +498,7 @@ TObjArray AliTPCTracklet::CreateTracklets(const AliTPCseed *s,
 
   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()];
   }
@@ -203,7 +507,188 @@ TObjArray AliTPCTracklet::CreateTracklets(const AliTPCseed *s,
   TObjArray tracklets;
   if (maxTracklets>72) maxTracklets=72; // just to protect against "users".
   for (Int_t i=0;i<maxTracklets&&sectors[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();
+    }
+  }
+  */
+}
index 4a89af7f195be9037f1a25d5a03f9e75155659e3..efaab0d32775de25a2d7f4e64c53a6fc6d1c9e59 100755 (executable)
 class TObjArray;
 class AliTPCseed;
 class AliExternalTrackParam;
+class AliTPCclusterMI;
+
+#include "TEllipse.h"
 
 class AliTPCTracklet:public TObject {
-public:
+public: 
+  enum TrackType {kKalman,kRiemann,kLinear,kQuadratic};
+
   AliTPCTracklet();
-  AliTPCTracklet(const AliTPCseed *s,Int_t sector);
+  AliTPCTracklet(const AliTPCseed *s,Int_t sector,TrackType type=kKalman,
+                Bool_t storeClusters=kFALSE);
+  AliTPCTracklet(const TObjArray &clusters,Int_t sector,TrackType type=kKalman,
+                Bool_t storeClusters=kFALSE);
   AliTPCTracklet(const AliTPCTracklet &t);
   AliTPCTracklet& operator=(const AliTPCTracklet &t);
   virtual ~AliTPCTracklet();
 
+  static TObjArray CreateTracklets(const TObjArray &clusters,
+                                  TrackType type=kKalman,
+                                  Bool_t storeClusters=kFALSE,
+                                  Int_t minClusters=0,
+                                  Int_t maxTracklets=72);
+
   static TObjArray CreateTracklets(const AliTPCseed *s,
+                                  TrackType type=kKalman,
+                                  Bool_t storeClusters=kFALSE,
                                   Int_t minClusters=0,
                                   Int_t maxTracklets=72);
 
+  static Bool_t PropagateToMeanX(const AliTPCTracklet &t1,
+                                const AliTPCTracklet &t2,
+                                AliExternalTrackParam *&t1m,
+                                AliExternalTrackParam *&t2m);
+
   // Returns the tracklet parametrisation at its outer most cluster.
   AliExternalTrackParam* GetOuter() const {return fOuter;};
   // Returns the tracklet parametrisation at its inner most cluster.
@@ -38,8 +59,33 @@ public:
   Int_t GetSector() const {return fSector;}
   // Returns the number of clusters assined to the tracklet.
   Int_t GetNClusters() const {return fNClusters;}
+  // Returns the clusters of this tracklet. In case they weren't stored it
+  // returns 0.
+  AliTPCclusterMI* GetClusters() const {return fClusters;};
+  // Test the functionality of the class. Generates some random tracks and
+  // refits them into tracklets. 
+  static void Test(const char *filename);
+  static void RandomND(Int_t ndim,const Double_t *p,const Double_t *c,
+                      Double_t *x);
+  static TEllipse ErrorEllipse(Double_t x,Double_t y,
+                              Double_t sx,Double_t sy,Double_t sxy);
 private:
+  static const Double_t kB2C; //! ugly to have the track parametrised in a way, that constand is allways needed
+  static double GetBz(Double_t *xyz);
+  void FitLinear(const AliTPCseed *track,Int_t sector,TrackType type);
+  void FitKalman(const AliTPCseed *track,Int_t sector);
+  void FitRiemann(const AliTPCseed *track,Int_t sector);
+  void Quadratic2Helix(Double_t *a,Double_t *ca,
+                      Double_t *b,Double_t *cb,
+                      Double_t x0,
+                      Double_t *p,Double_t *c);
+  Bool_t Riemann2Helix(Double_t *a,Double_t *ca,
+                      Double_t *b,Double_t *cb,
+                      Double_t x0,
+                      Double_t *p,Double_t *c);   
   Int_t fNClusters; // The number of clusters assined to the tracklet.
+  Int_t fNStoredClusters; // The number of stored clusters.
+  AliTPCclusterMI *fClusters; //[fNStoredClusters] The clusters of the track, if stored (otherwise 0)
   Int_t fSector; // The sector this tracklet lives in.
   AliExternalTrackParam *fOuter; // The tracklet parametrisation at its outer most cluster.
   AliExternalTrackParam *fInner; // The tracklet parametrisation at its inner most cluster.
diff --git a/TPC/AliTPCcalibAlign.cxx b/TPC/AliTPCcalibAlign.cxx
new file mode 100644 (file)
index 0000000..8b1e8f2
--- /dev/null
@@ -0,0 +1,443 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//     Class to make a internal alignemnt of TPC chambers                    //
+//
+//     Different linear tranformation investigated
+//     12 parameters - arbitrary transformation 
+//      9 parameters - scaling fixed to 1
+//      6 parameters - 
+////
+//// Only the 12 parameter version works so far...
+////
+////
+//// 
+
+#include "TLinearFitter.h"
+#include "AliTPCcalibAlign.h"
+#include "AliExternalTrackParam.h"
+
+#include <iostream>
+using namespace std;
+
+ClassImp(AliTPCcalibAlign)
+
+AliTPCcalibAlign::AliTPCcalibAlign()
+  :fFitterArray12(72*72),fFitterArray9(72*72),fFitterArray6(72*72)
+{
+  //
+  // Constructor
+  //
+  for (Int_t i=0;i<72*72;++i) {
+    fPoints[i]=0;
+  }
+}
+
+AliTPCcalibAlign::~AliTPCcalibAlign() {
+  //
+  // destructor
+  //
+}
+
+void AliTPCcalibAlign::Process(const AliExternalTrackParam &tp1,
+                              const AliExternalTrackParam &tp2,
+                              Int_t s1,Int_t s2) {
+  //
+  // Process function to fill fitters
+  //
+  Double_t t1[5],t2[5];
+  Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
+  Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
+  x1   =tp1.GetX();
+  y1   =tp1.GetY();
+  z1   =tp1.GetZ();
+  Double_t snp1=tp1.GetSnp();
+  dydx1=snp1/TMath::Sqrt(1.-snp1*snp1);
+  Double_t tgl1=tp1.GetTgl();
+  // dz/dx = 1/(cos(theta)*cos(phi))
+  dzdx1=1./TMath::Sqrt((1.+tgl1*tgl1)*(1.-snp1*snp1));
+  //  x2  =tp1->GetX();
+  y2   =tp2.GetY();
+  z2   =tp2.GetZ();
+  Double_t snp2=tp2.GetSnp();
+  dydx2=snp2/TMath::Sqrt(1.-snp2*snp2);
+  Double_t tgl2=tp2.GetTgl();
+  dzdx1=1./TMath::Sqrt((1.+tgl2*tgl2)*(1.-snp2*snp2));
+
+  Process12(t1,t2,GetOrMakeFitter12(s1,s2));
+  Process9(t1,t2,GetOrMakeFitter12(s1,s2));
+  Process6(t1,t2,GetOrMakeFitter12(s1,s2));
+  ++fPoints[72*s1+s2];
+}
+
+void AliTPCcalibAlign::Process12(Double_t *t1,
+                                Double_t *t2,
+                                TLinearFitter *fitter) {
+  // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
+  // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
+  // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
+  // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
+  // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
+  //
+  //       a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
+  //       a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
+  //       a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11] 
+  Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
+  Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
+
+  // TODO:
+  Double_t sy    = 1.;
+  Double_t sz    = 1.;
+  Double_t sdydx = 1.;
+  Double_t sdzdx = 1.;
+
+  Double_t p[12];
+  Double_t value;
+
+  // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
+  // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
+  // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
+  for (Int_t i=0; i<12;i++) p[i]=0.;
+  p[3+0] = x1;          // a10
+  p[3+1] = y1;          // a11
+  p[3+2] = z1;          // a12
+  p[9+1] = 1.;          // a13
+  p[0+1] = y1*dydx2;    // a01
+  p[0+2] = z1*dydx2;    // a02
+  p[9+0] = dydx2;       // a03
+  value  = y2;
+  fitter->AddPoint(p,value,sy);
+
+  // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
+  // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
+  // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 + a03)*dzdx2;
+  for (Int_t i=0; i<12;i++) p[i]=0.;
+  p[6+0] = x1;           // a20 
+  p[6+1] = y1;           // a21
+  p[6+2] = z1;           // a22
+  p[9+2] = 1.;           // a23
+  p[0+1] = y1*dzdx2;     // a01
+  p[0+2] = z1*dzdx2;     // a02
+  p[9+0] = dzdx2;        // a03
+  value  = z2;
+  fitter->AddPoint(p,value,sz);
+
+  // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
+  // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
+  for (Int_t i=0; i<12;i++) p[i]=0.;
+  p[3+0] = 1.;           // a10
+  p[3+1] = dydx1;        // a11
+  p[3+2] = dzdx1;        // a12
+  p[0+0] = -dydx2;       // a00
+  p[0+1] = -dydx1*dydx2; // a01
+  p[0+2] = -dzdx1*dydx2; // a02
+  value  = 0.;
+  fitter->AddPoint(p,value,sdydx);
+
+  // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
+  // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
+  for (Int_t i=0; i<12;i++) p[i]=0.;
+  p[6+0] = 1;            // a20
+  p[6+1] = dydx1;        // a21
+  p[6+2] = dzdx1;        // a22
+  p[0+0] = -dzdx2;       // a00
+  p[0+1] = -dydx1*dzdx2; // a01
+  p[0+2] = -dzdx1*dzdx2; // a02
+  value  = 0.;
+  fitter->AddPoint(p,value,sdzdx);
+}
+
+void AliTPCcalibAlign::Process9(Double_t *t1,
+                               Double_t *t2,
+                               TLinearFitter *fitter) {
+  // x2    =  a00*x1 + a01*y1 + a02*z1 + a03
+  // y2    =  a10*x1 + a11*y1 + a12*z1 + a13
+  // z2    =  a20*x1 + a21*y1 + a22*z1 + a23
+  // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
+  // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
+  //
+  //       a00  a01 a02  a03     p[0]   p[1]  p[2]  p[9]
+  //       a10  a11 a12  a13 ==> p[3]   p[4]  p[5]  p[10]
+  //       a20  a21 a22  a23     p[6]   p[7]  p[8]  p[11] 
+  Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
+  Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
+
+  // TODO:
+  Double_t sy    = 1.;
+  Double_t sz    = 1.;
+  Double_t sdydx = 1.;
+  Double_t sdzdx = 1.;
+
+  Double_t p[12];
+  Double_t value;
+
+  // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
+  // y2  =  a10*x1 + a11*y1 + a12*z1 + a13
+  // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a03)*dydx2
+  for (Int_t i=0; i<12;i++) p[i]=0.;
+  p[3+0] = x1;
+  p[3+1] = y1;
+  p[3+2] = z1;
+  p[9+1] = 1.;
+  p[0+1] = y1*dydx2;
+  p[0+2] = z1*dydx2;
+  p[9+0] = dydx2;
+  value  = y2;
+  fitter->AddPoint(p,value,sy);
+
+  // x2  =  a00*x1 + a01*y1 + a02*z1 + a03
+  // z2  =  a20*x1 + a21*y1 + a22*z1 + a23
+  // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a03)*dzdx2;
+  for (Int_t i=0; i<12;i++) p[i]=0.;
+  p[6+0] = x1;
+  p[6+1] = y1;
+  p[6+2] = z1;
+  p[9+2] = 1.;
+  p[0+1] = y1*dzdx2;
+  p[0+2] = z1*dzdx2;
+  p[9+0] = dzdx2;
+  value  = z2;
+  fitter->AddPoint(p,value,sz);
+
+  // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
+  // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
+  for (Int_t i=0; i<12;i++) p[i]=0.;
+  p[3+0] = 1.;
+  p[3+1] = dydx1;
+  p[3+2] = dzdx1;
+  p[0+0] = -dydx2;
+  p[0+1] = -dydx1*dydx2;
+  p[0+2] = -dzdx1*dydx2;
+  value  = 0.;
+  fitter->AddPoint(p,value,sdydx);
+
+  // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
+  // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
+  for (Int_t i=0; i<12;i++) p[i]=0.;
+  p[6+0] = 1;
+  p[6+1] = dydx1;
+  p[6+2] = dzdx1;
+  p[0+0] = -dzdx2;
+  p[0+1] = -dydx1*dzdx2;
+  p[0+2] = -dzdx1*dzdx2;
+  value  = 0.;
+  fitter->AddPoint(p,value,sdzdx);
+}
+
+void AliTPCcalibAlign::Process6(Double_t *t1,
+                               Double_t *t2,
+                               TLinearFitter *fitter) {
+  // x2    =  1  *x1 +-a01*y1 + 0      +a03
+  // y2    =  a01*x1 + 1  *y1 + 0      +a13
+  // z2    =  a20*x1 + a21*y1 + 1  *z1 +a23
+  // dydx2 = (a10    + a11*dydx1 + a12*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
+  // dzdx2 = (a20    + a21*dydx1 + a22*dzdx1)/(a00    + a01*dydx1 + a02*dzdx1)
+  //
+  //       1   -a01 0    a03     x     -p[0]  x     p[3]
+  //       a10  1   0    a13 ==> p[0]   x     x     p[4]
+  //       a20  a21 1    a23     p[1]   p[2]  x     p[5] 
+  Double_t &x1=t1[0], &y1=t1[1], &z1=t1[2], &dydx1=t1[3], &dzdx1=t1[4];
+  Double_t &x2=t2[0], &y2=t2[1], &z2=t2[2], &dydx2=t2[3], &dzdx2=t2[4];
+
+  // TODO:
+  Double_t sy    = 1.;
+  Double_t sz    = 1.;
+  Double_t sdydx = 1.;
+  Double_t sdzdx = 1.;
+
+  Double_t p[12];
+  Double_t value;
+
+  // x2  =  1  *x1 +-a01*y1 + 0      +a03
+  // y2  =  a01*x1 + 1  *y1 + 0      +a13
+  // y2' =  a10*x1 + a11*y1 + a12*z1 + a13 + (a01*y1 + a02*z1 + a03)*dydx2
+  for (Int_t i=0; i<12;i++) p[i]=0.;
+  p[3+0] = x1;
+  p[3+1] = y1;
+  p[3+2] = z1;
+  p[9+1] = 1.;
+  p[0+1] = y1*dydx2;
+  p[0+2] = z1*dydx2;
+  p[9+0] = dydx2;
+  value  = y2;
+  fitter->AddPoint(p,value,sy);
+
+  // x2  =  1  *x1 +-a01*y1 + 0      + a03
+  // z2  =  a20*x1 + a21*y1 + 1  *z1 + a23
+  // z2' =  a20*x1 + a21*y1 + a22*z1 + a23 + (a01*y1 + a02*z1 +a03)*dzdx2;
+  for (Int_t i=0; i<12;i++) p[i]=0.;
+  p[6+0] = x1;
+  p[6+1] = y1;
+  p[6+2] = z1;
+  p[9+2] = 1.;
+  p[0+1] = y1*dzdx2;
+  p[0+2] = z1*dzdx2;
+  p[9+0] = dzdx2;
+  value  = z2;
+  fitter->AddPoint(p,value,sz);
+
+  // dydx2 = (a10 + a11*dydx1 + a12*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
+  // (a10 + a11*dydx1 + a12*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dydx2 = 0
+  for (Int_t i=0; i<12;i++) p[i]=0.;
+  p[3+0] = 1.;
+  p[3+1] = dydx1;
+  p[3+2] = dzdx1;
+  p[0+0] = -dydx2;
+  p[0+1] = -dydx1*dydx2;
+  p[0+2] = -dzdx1*dydx2;
+  value  = 0.;
+  fitter->AddPoint(p,value,sdydx);
+
+  // dzdx2 = (a20 + a21*dydx1 + a22*dzdx1)/( a00 + a01*dydx1 + a02*dzdx1)
+  // (a20 + a21*dydx1 + a22*dzdx1) - (a00 + a01*dydx1 + a02*dzdx1)*dzdx2 = 0
+  for (Int_t i=0; i<12;i++) p[i]=0.;
+  p[6+0] = 1;
+  p[6+1] = dydx1;
+  p[6+2] = dzdx1;
+  p[0+0] = -dzdx2;
+  p[0+1] = -dydx1*dzdx2;
+  p[0+2] = -dzdx1*dzdx2;
+  value  = 0.;
+  fitter->AddPoint(p,value,sdzdx);
+}
+
+void AliTPCcalibAlign::Eval() {
+  TLinearFitter *f;
+  for (Int_t s1=0;s1<72;++s1)
+    for (Int_t s2=0;s2<72;++s2)
+      if ((f=GetFitter12(s1,s2))&&fPoints[72*s1+s2]>12) {
+       //      cerr<<s1<<","<<s2<<": "<<fPoints[72*s1+s2]<<endl;
+       if (!f->Eval()) {
+         cerr<<"Evaluation failed for "<<s1<<","<<s2<<endl;
+       }
+      }
+  /*
+                   
+  fitter->Eval();
+  fitter->Eval();
+  chi212 = align->GetChisquare()/(4.*entries);
+
+  TMatrixD mat(13,13);
+  TVectorD par(13);
+  align->GetParameters(par);
+  align->GetCovarianceMatrix(mat);
+
+  //
+  //
+  for (Int_t i=0; i<12;i++){
+    palign12(i)= par(i+1);
+    for (Int_t j=0; j<12;j++){
+      pcovar12(i,j)   = mat(i+1,j+1);
+      pcovar12(i,j) *= chi212;
+    }
+  }
+  //
+  for (Int_t i=0; i<12;i++){
+    psigma12(i)  = TMath::Sqrt(pcovar12(i,i));
+    palignR12(i) = palign12(i)/TMath::Sqrt(pcovar12(i,i));
+    for (Int_t j=0; j<12;j++){
+      pcovarN12(i,j) = pcovar12(i,j)/TMath::Sqrt(pcovar12(i,i)*pcovar12(j,j));
+    }
+  }
+  */
+}
+
+Bool_t AliTPCcalibAlign::GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a) {
+  if (!GetFitter12(s1,s2))
+    return false;
+  else {
+    TVectorD p(12);
+    cerr<<"foo"<<endl;
+    GetFitter12(s1,s2)->GetParameters(p);
+    cerr<<"bar"<<endl;
+    a.ResizeTo(4,4);
+    a[0][0]=p[0];
+    a[0][1]=p[1];
+    a[0][2]=p[2];
+    a[0][3]=p[9];
+    a[1][0]=p[3];
+    a[1][1]=p[4];
+    a[1][2]=p[5];
+    a[1][3]=p[10];
+    a[2][0]=p[6];
+    a[2][1]=p[7];
+    a[2][2]=p[8];
+    a[2][3]=p[11];
+    a[3][0]=0.;
+    a[3][1]=0.;
+    a[3][2]=0.;
+    a[3][3]=1.;
+    return true;
+  } 
+}
+
+Bool_t AliTPCcalibAlign::GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a) {
+  if (!GetFitter9(s1,s2))
+    return false;
+  else {
+    TVectorD p(9);
+    GetFitter9(s1,s2)->GetParameters(p);
+    a.ResizeTo(4,4);
+    a[0][0]=p[0];
+    a[0][1]=p[1];
+    a[0][2]=p[2];
+    a[0][3]=p[9];
+    a[1][0]=p[3];
+    a[1][1]=p[4];
+    a[1][2]=p[5];
+    a[1][3]=p[10];
+    a[2][0]=p[6];
+    a[2][1]=p[7];
+    a[2][2]=p[8];
+    a[2][3]=p[11];
+    a[3][0]=0.;
+    a[3][1]=0.;
+    a[3][2]=0.;
+    a[3][3]=1.;
+    return true;
+  } 
+}
+
+Bool_t AliTPCcalibAlign::GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a) {
+  if (!GetFitter6(s1,s2))
+    return false;
+  else {
+    TVectorD p(6);
+    cerr<<"foo"<<endl;
+    GetFitter6(s1,s2)->GetParameters(p);
+    cerr<<"bar"<<endl;
+    a.ResizeTo(4,4);
+    a[0][0]=p[0];
+    a[0][1]=p[1];
+    a[0][2]=p[2];
+    a[0][3]=p[9];
+    a[1][0]=p[3];
+    a[1][1]=p[4];
+    a[1][2]=p[5];
+    a[1][3]=p[10];
+    a[2][0]=p[6];
+    a[2][1]=p[7];
+    a[2][2]=p[8];
+    a[2][3]=p[11];
+    a[3][0]=0.;
+    a[3][1]=0.;
+    a[3][2]=0.;
+    a[3][3]=1.;
+    return true;
+  } 
+}
diff --git a/TPC/AliTPCcalibAlign.h b/TPC/AliTPCcalibAlign.h
new file mode 100644 (file)
index 0000000..192e376
--- /dev/null
@@ -0,0 +1,75 @@
+#ifndef ALITPCCALIBALIGN_H
+#define ALITPCCALIBALIGN_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+////
+////
+////
+
+#include "TObject.h"
+#include "TObjArray.h"
+#include "TLinearFitter.h"
+
+class AliExternalTrackParam;
+
+class AliTPCcalibAlign:public TObject {
+public:
+  AliTPCcalibAlign();
+
+  virtual ~AliTPCcalibAlign();
+
+  void Process(const AliExternalTrackParam &t1,
+              const AliExternalTrackParam &t2,
+              Int_t s1,Int_t s2);
+  void Eval();
+  TLinearFitter* GetFitter12(Int_t s1,Int_t s2) {
+    return static_cast<TLinearFitter*>(fFitterArray12[s1*72+s2]);
+  }
+  TLinearFitter* GetFitter9(Int_t s1,Int_t s2) {
+    return static_cast<TLinearFitter*>(fFitterArray9[s1*72+s2]);
+  }
+  TLinearFitter* GetFitter6(Int_t s1,Int_t s2) {
+    return static_cast<TLinearFitter*>(fFitterArray6[s1*72+s2]);
+  }
+  Bool_t GetTransformation12(Int_t s1,Int_t s2,TMatrixD &a);
+  Bool_t GetTransformation9(Int_t s1,Int_t s2,TMatrixD &a);
+  Bool_t GetTransformation6(Int_t s1,Int_t s2,TMatrixD &a);
+private:
+  void Process12(Double_t *t1,
+                Double_t *t2,
+                TLinearFitter *fitter);
+  void Process9(Double_t *t1,
+               Double_t *t2,
+               TLinearFitter *fitter);
+  void Process6(Double_t *t1,
+               Double_t *t2,
+               TLinearFitter *fitter);
+  TLinearFitter* GetOrMakeFitter12(Int_t s1,Int_t s2) {
+    //get or make fitter
+    if (!fFitterArray12[s1*72+s2])
+      fFitterArray12[s1*72+s2]=new TLinearFitter(12,"x0++x1++x2++x3++x4++x5++x6++x7++x8++x9++x10++x11");
+    return GetFitter12(s1,s2);
+  }
+  TLinearFitter* GetOrMakeFitter9(Int_t s1,Int_t s2) {
+    //get or make fitter
+    if (!fFitterArray9[s1*72+s2])
+      fFitterArray9[s1*72+s2]=new TLinearFitter(9,"x0++x1++x2++x3++x4++x5++x6++x7++x8");
+    return GetFitter9(s1,s2);
+  }
+  TLinearFitter* GetOrMakeFitter6(Int_t s1,Int_t s2) {
+    //get or make fitter
+    if (!fFitterArray6[s1*72+s2])
+      fFitterArray6[s1*72+s2]=new TLinearFitter(6,"x0++x1++x2++x3++x4++x5");
+    return GetFitter6(s1,s2);
+  }
+  TObjArray fFitterArray12;
+  TObjArray fFitterArray9;
+  TObjArray fFitterArray6;
+  Int_t fPoints[72*72];
+
+  ClassDef(AliTPCcalibAlign,1)
+};
+
+#endif
diff --git a/TPC/AliTPCcalibAlignment.cxx b/TPC/AliTPCcalibAlignment.cxx
new file mode 100644 (file)
index 0000000..dbe2e40
--- /dev/null
@@ -0,0 +1,123 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//     Class to make and monitor alignment of the TPC                        //
+// 
+//         
+
+
+#include "AliTracker.h"
+#include "AliTPCcalibAlignment.h"
+#include "AliTPCseed.h"
+#include "AliTPCclusterMI.h"
+#include "AliTPCClusterParam.h"
+#include "AliESDtrack.h"
+#include "TTreeStream.h"
+#include "AliTPCTracklet.h"
+
+#include <iostream>
+using namespace std;
+
+/* TEST
+
+  gSystem->Load("$ALICE_ROOT/TPC/TPCcalib/libTPCcalib.so");
+
+  .L AliXRDPROOFtoolkit.cxx+
+  AliXRDPROOFtoolkit tool;
+  TChain * chain = tool.MakeChain("listcosmic.txt","esdTree",0,100,0)
+//TChain * chain = tool.MakeChain("listpp.txt","esdTree",0,100,0)
+
+  AliTPCcalibTracks::AddInfo(chain,"TPCClusterParam.root");
+  AliTPCcalibTracks::AddCuts(chain,"LOWFLUX");
+
+//  AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 0., 10., 2); // NO FIELD
+  AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., 2); // FIELD
+//  field->SetL3ConstField(0); //Using const. field in the barrel
+  AliTracker::SetFieldMap(field,kFALSE);
+
+   chain->SetBranchStatus("*",1);
+  chain->Process("$ALICE_ROOT/TPC/TPCcalib/AliTPCSelectorTracks.cxx+"); 
+
+  TFile f("AliTPCcalibAlignmentDebug.root")  
+  TTree *t=f.Get("Alignment")
+
+  t->Draw("tracklet1.fInner.GetY()-trackletR1.fInner.GetY()")
+  t->Draw("tracklet1.fInner.GetSnp()-tracklet1.fOuter.GetSnp()")
+  t->Draw("common1.GetSnp()-common2.GetSnp()")
+  t->Draw("common1.GetY()-common2.GetY()")
+  t->Draw("common1.GetZ()-common2.GetZ()")
+
+  t->Draw("track1.GetY()-track2.GetY()","abs(track1.GetY()-track2.GetY())<0.2")
+  t->Draw("track1.GetZ()-track2.GetZ()","abs(track1.GetZ()-track2.GetZ())<0.2")
+  t->Draw("track1.GetY()-track2.GetY():track1.GetZ()-track2.GetZ()","abs(track1.GetY()-track2.GetY())<0.2&&abs(track1.GetZ()-track2.GetZ())<0.2")
+
+
+ */
+
+ClassImp(AliTPCcalibAlignment)
+
+AliTPCcalibAlignment::AliTPCcalibAlignment() {
+  //
+  // Constructor
+  //
+   fDebugStream=new TTreeSRedirector("AliTPCcalibAlignmentDebug.root");
+}
+
+AliTPCcalibAlignment::~AliTPCcalibAlignment() {
+  //
+  // Destructor
+  //
+  delete fDebugStream;
+}
+
+void AliTPCcalibAlignment::Process(AliTPCseed *track) {
+  //
+  // Track processing
+  //
+  TObjArray tracklets=
+    AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kKalman,kFALSE,20,2);
+  TObjArray trackletsL=
+    AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kLinear,kFALSE,20,2);
+  TObjArray trackletsQ=
+    AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kQuadratic,kFALSE,20,2);
+  TObjArray trackletsR=
+    AliTPCTracklet::CreateTracklets(track,AliTPCTracklet::kRiemann,kFALSE,20,2);
+  tracklets.SetOwner();
+  if (tracklets.GetEntries()==2) {
+    AliTPCTracklet* t1=static_cast<AliTPCTracklet*>(tracklets[0]);
+    AliTPCTracklet* t2=static_cast<AliTPCTracklet*>(tracklets[1]);
+    AliExternalTrackParam *common1=0,*common2=0;
+    if (AliTPCTracklet::PropagateToMeanX(*t1,*t2,common1,common2))
+      (*fDebugStream)<<"Alignment"
+                    <<"common1.="<<common1
+                    <<"common2.="<<common2
+                    <<"tracklet1.="<<tracklets[0]
+                    <<"tracklet2.="<<tracklets[1]
+                    <<"trackletL1.="<<trackletsL[0]
+                    <<"trackletL2.="<<trackletsL[1]
+                    <<"trackletQ1.="<<trackletsQ[0]
+                    <<"trackletQ2.="<<trackletsQ[1]
+                    <<"trackletR1.="<<trackletsR[0]
+                    <<"trackletR2.="<<trackletsR[1]
+                    <<"\n";
+    delete common1;
+    delete common2;
+ }
+}
+
+
diff --git a/TPC/AliTPCcalibAlignment.h b/TPC/AliTPCcalibAlignment.h
new file mode 100644 (file)
index 0000000..398b7a7
--- /dev/null
@@ -0,0 +1,29 @@
+#ifndef ALI_TPC_CALIB_ALIGNMENT_H
+#define ALI_TPC_CALIB_ALIGNMENT_H
+
+#include "TNamed.h"
+// ugly, but...
+class AliTPCseed;
+class AliExternalTrackParam;
+class AliESDtrack;
+class AliRieman;
+class AliTPCtrackerMI;
+class TTreeSRedirector;
+
+class AliTPCcalibAlignment:public TNamed {
+public:
+  AliTPCcalibAlignment();
+  virtual ~AliTPCcalibAlignment();
+
+  virtual void Process(AliTPCseed *track);
+
+private:
+  AliTPCcalibAlignment(const AliTPCcalibAlignment&);
+  AliTPCcalibAlignment& operator=(const AliTPCcalibAlignment&);
+
+  TTreeSRedirector *fDebugStream;
+
+  ClassDef(AliTPCcalibAlignment,1);
+};
+
+#endif
index 589ad415379ed4309d6d8236d96ce16fcd507164..cbe26e523eae579349e008317beb2315fc94c3c4 100644 (file)
@@ -20,7 +20,8 @@
 #pragma link C++ class  AliTPCSelectorTracks+;
 #pragma link C++ class  AliAnaTPCTrackBase+;
 #pragma link C++ class  AliAnaTPCTrackCalib+;
-
+#pragma link C++ class  AliTPCcalibAlign+;
+#pragma link C++ class  AliTPCcalibAlignment+;
 
 #endif
 
index 3a423cf75d8b687100dc93a08873ba4154c00005..7cb3dc4330883e8175be3057c159ade8c0ecd521 100644 (file)
@@ -1,6 +1,6 @@
 
 SRCS = AliTPCcalibTracksCuts.cxx   AliTPCcalibTracks.cxx   AliTPCcalibTracksGain.cxx  \
-        AliTPCSelectorESD.cxx   AliTPCSelectorTracks.cxx   AliTPCCalPadRegion.cxx AliTPCFitPad.cxx AliAnaTPCTrackBase.cxx AliAnaTPCTrackCalib.cxx
+        AliTPCSelectorESD.cxx   AliTPCSelectorTracks.cxx   AliTPCCalPadRegion.cxx AliTPCFitPad.cxx AliAnaTPCTrackBase.cxx AliAnaTPCTrackCalib.cxx AliTPCcalibAlign.cxx AliTPCcalibAlignment.cxx 
       
 
 HDRS:= $(SRCS:.cxx=.h)